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1. Introduction 



The idea that extra spatial dimensions may provide a solution to the hierarchy problem 

by lowering the scale of gravity to the TcV region [1-5] also raises the exciting possibility 
of black hole (BH) production in elementary particle collisions at high energies [6,7]. For 
a sufficient number of "large" extra dimensions (three or more), the current experimental 
and astrophysical limits [8,9] do not rule out such a process occurring at the energy scales 
accessible to the Large Hadron Collider (LHC). Even if not realized at the LHC, the 
production and decay of extra-dimensional black holes poses an important problem in 
theoretical physics. Given the number of dimensions D = 4+n and the fundamental Planck 
scale Md, one should be able to describe the production of black holes at collision energies 
well above Md using classical general relativity extended to (4 + n)-dimensions [10-13]. In 
addition, at least at black hole masses well above Md, the decay of a black hole so formed 
should be describable in terms of the Hawking radiation [14, 15] expected from such an 
object. 

Although many important questions remain unanswered, there has been steady progress 
in the theoretical understanding of these issues in recent years. Limits have been placed 
on the fractions of the collision energy and angular momentum that can end up in the 
black hole, as functions of the impact parameter of the collision [18-22]. The fact that 
the black hole will in general have non-zero angular momentum has led to detailed stud- 
ies of the Hawking radiation from rotating higher-dimensional black holes [23-34], which 
have revealed interesting features such as differences in the angular distributions of emitted 
particles with different spins, as well as modifications of the energy distributions. 

In parallel with the purely theoretical work, efforts have been made to incorporate the 
results obtained into Monte Carlo programs [7, 35-38] that generate simulated black hole 
events, for use in studies of the experimental detectability of this process. In the present 
paper we report on progress in refining the widely-used CHARYBDIS event generator [35, 
39,40], to take into account the recent theoretical work on the production and decay of 
rotating black holes. This enables us to study the effects of rotation on experimentally 
observable quantities, assuming that black hole production does indeed take place. 

In the following section we summarize the current theoretical results on the formation 
of black holes in particle collisions and the implied constraints on their masses and angular 
momenta. We also describe how these results are incorporated into the CHARYBDIS2^ 
simulation, with various models for the joint probability distribution of mass and spin 
within the allowed region. Then, in Sect. 3, we discuss the spin-down and decay of the 
black hole through Hawking emission of Standard-Model (SM) particles confined to the 
physical 3-brane. Here again the simulation is refined to take full account of the theoretical 
results, and options are included for the modelling of aspects that are not well understood, 
such as back-reaction effects. 

In any simulation of black hole decay in theories with low-scale gravity, one eventually 
reaches the stage at which the mass and/or temperature of the black hole are comparable 

^We will refer from now on to the new release as CHARYBDIS2 and will reserve CHARYBDIS for earlier 
versions. The particular version described here is CHARYBDIS2 . 0. 
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with the Planck scale, whereupon a complete theory of quantum gravity is required. In the 
absence of such a theory, various models for the termination of black hole decay have been 
suggested [41-46]. In Sect. 4 we explain the extended range of options we have included 
for this phase of the simulation. 

Section 5 presents some results of the new CHARYBDIS2 simulation, with particular 
emphasis on consequences of black hole rotation. Our purpose here is not to perform a 
detailed experimental analysis, but to illustrate new features that will need to be taken into 
account in attempts to extract fundamental parameters from experimental data, should this 
process actually occur. 

Finally in Sect. 6 we draw some conclusions and make suggestions for further work. 
We also discuss the comparative features of CHARYBDIS2 and the black hole event generator 
BlackMax [37,38], which is currently the only other simulation that takes black hole rotation 
into account. 



2. Black hole production 

Production of a black hole in a hadron-hadron collision occurs when two of the colliding 

partons pass sufficiently closely that they become trapped by their mutual gravitational 
attraction. Prior to the collision, we may take these partons to be travelling in anti-parallel 
directions, with their directions of motion separated by an impact parameter b. A complex 
shaped event horizon forms, which must quickly relax to one of the axisymmetric stationary 
black hole solutions of Einstein's equations according to the 'no hair' theorem of classical 
general relativity [15]. For that reason this is sometimes called the balding phase. An 
emission of gravitational and gauge radiation accompanies the production event. 

It has been argued [6] that the typical time scale for loss of such asymmetries is 
of the order of the horizon radius rn (note we are using natural units - see appendix 
A). This is physically reasonable because the asymmetries are related to a distortion 
of the geometry (with respect to the stationary solution); regardless of the details of the 
interaction responsible for removing them, it will involve a signal propagating over a region 
of typical size rn- On the other hand the time scale At for evaporation of a black hole with 
a mass well above the Planck scale is controlled by the Hawking energy flux. Since we are 
only interested in an order of magnitude, we can estimate it by using the D-dimensional 
Schwarzschild case combined with dimensional analysis (see for example Sect. 3 of [6]) to 
obtain 

At f M \ ^ , , 

0^ ■ (2-1) 



The prefactor is a constant of order unity containing the dimensionally reduced energy flux 
and other convention-dependent constants, and Md is the Planck mass in D dimensions^. 
For M 3> Md the time for evaporation will be much longer than for balding. 

For the theoretical calculations involved in the model for formation and for the Hawking 
fluxes used in the evaporation, we assume that the gauge charges of the incoming partons 



^See Appendix A for our convention on the definition of Mb- 



-3- 



are completely discharged during the production/balding phase by Schwinger emission [16, 
17]. Further, we assume that the black hole solution formed is always of the Myers-Perry 
type [47] (rather than the alternative 'black ring' type, which is known to exist in five 
dimensions [48], and for which there is strong evidence for D > 6 [49]). Given these 
assumptions, two questions are relevant to the theoretical modelling of the production 
phase. The first is how small the impact parameter b has to be before two partons will 
produce a black hole, or equivalently, what is the parton-level cross section for black hole 
formation. The second is what fractions of the initial state mass and angular momentum 
are trapped within the Myers-Perry black hole during the production phase. 

Theoretical techniques used in recent years have provided more rigorous and complete 
answers to these questions than were available at the last release of CHARYBDIS. Some 
of these techniques are discussed in Sect. 2.1. Their incorporation into the program, to 
produce a more accurate simulation of the production phase, is discussed in Sect. 2.2. 

2.1 Theoretical studies of the production phase 

In the production phase, the system may be described reasonably well using classical 
physics, provided the parton collision energy is sufficiently far above the Planck scale [10- 
13]. In principle, one could find reasonably accurate answers to the questions posed above 
by solving for the spacctimc in the fiiturc of two colliding partons. However, this has not 
been achieved, even numerically in simplified scenarios. 

The standard approach to studying the production phase has been the trapped surface 
method [18-22]. This method utilises the fact that the black hole horizon begins forming in 
the spacetime region outside the future lightcone of the collision event, where we can solve 
for the geometry. By studying apparent horizons in spacetime slices of this region, and in 
particular looking at certain areas associated with the apparent horizon as b is varied, one 
can set bounds on the parton-level cross section and the mass M and angular momentum 
J trapped for given b. 

Few trapped surface calculations conducted thus far have attempted to obtain results 
for nonzero b, making them of limited use to our simulation. The most detailed that has, 
is that of Yoshino and Rychkov [22] who considered all of the numbers of dimensions used 
by CHARYBDIS2. They produced (M, J) bounds up to the maximum impact parameters 
which give rise to apparent horizons in their method. The results of this calculation are 
therefore the primary theoretical input to our simulation of the production phase. A more 
detailed discussion of the Yoshino-Rychkov method is given in Appendix B. 

It should be noted that the Yoshino-Rychkov calculation models the colliding partons 
as boosted Schwarzschild-Tangherlini black holes, and so neglects the effects of the spin, 
charge, and finite size of the colliding partons. For each effect individually, trapped surface 
calculations have been carried out [50-52] - but all are only for 6 = 0. One important 
observation from the calculations is that charge effects may be significant. These will be 
included in CHARYBDIS2 when the relevant calculations are extended to nonzero impact 
parameter. 

Alongside the trapped surface method, alternative techniques have been developed 
which use a perturbative approach and/or other approximations to estimate directly the 
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Figure 1: A comparison of various theoretical results for the production phase mass loss in the 
6 = case, with the average mass loss produced for 6 = by our simulation. Black squares: 
'Particle falling into black hole' results from [53]. Open squares: 'Instantaneous collision' results 
from [54]. Asterisk: 'Bondi's news function' result from [55-57]. Crosses: 'Trapped surface method' 
upper bound on mass loss. Points with error bars: average 6 = mass loss from our simulation. 
The error bars represent the standard deviation in the 6 = output. 



mass lost in the production phase. In one setup [53], the collision is modelled as an 
ultra-relativistic particle falling into a Schwarzschild-Tangherlini black hole, and the grav- 
itational emission is calculated by assuming that the gravitational effects of the in-falling 
particle may be treated as a perturbation on top of the Schwarzschild-Tangherlini metric. 
Another calculation [54] also uses a perturbative approach and assumes that the collision 
is instantaneous. As a final example, D'Eath and Payne [55-57] have estimated the mass 
loss in the D = A axisymmetric collision case by finding the first two terms of Bondi's 
news function, and then extrapolating off axis. Here some assumptions about the angular 
dependence of the radiation are made. 

The results produced so far by these methods have been limited to 6 = and certain 
values of D. The 6 = results from different techniques arc compared with the trapped 
surface bound in Fig. 1. The general indication from these is that much less mass is lost 
during the production phase than the Yoshino-Rychkov upper limits indicate. 

2.2 Incorporation of the results into CHARYBDIS2 
Cross section 

In earlier versions of CHARYBDIS, parton-level cross sections for different D values were 
calculated according to the simple formula a = Trrg{-\/s) which is based on Thome's hoop 
conjecture [58]. Here rs{^/s) is the radius of the Z)-dimensional Schwarzschild black hole 
with mass ^/s. Incorporation of the Yoshino-Rychkov cross section results simply requires 
multiplying these a values by the 'formation factors' given in Table II of [22]. The increase 
in a ranges from a factor ofl.5atD = 5to3.2atD = 11. The maximum impact 
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parameters for black hole production, bmaxi to be generated in CHARYBDIS2, is adjusted 
accordingly (the two are related through a = tt^^^^.). 

Mass and angular momentum loss 

Following Yoshino and Rychkov, we denote the fractions of the initial state mass and 
angular momentum, trapped after production, by ^ and ^ respectively. For a given number 
of dimensions D and impact parameter b, the Yoshino-Rychkov bound on these quantities 
is a curve in the (^, (") plane. Examples of such curves for various D and b values are 
given as the solid lines in Fig. 2. A boundary curve S,h{C) always possesses the following 
two key properties. First, it always passes through C = with a ^ value between and 
1, where ^(,(0) = ^n,{b,D) in the language of [22]. Second, it increases monotonically after 
this, passing through ( = 1 with a value satisfying ^u, < < 1- The allowed region is 
then delimited by this curve and the lines C = 0, ^ = 1, and ^ = 1. 

The new simulation of mass and angular momentum loss based on these curves consists 
of a point being generated at random in the square 0<.^<1,0<C^1- The probabil- 
ity distribution for generating this point goes to zero along the Yoshino-Rychkov bound 
corresponding to the D and 6 values of the event, such that the generated point is always 
inside the bound. The ^ and ( coordinates of this point are then taken as the fractional 
mass and angular momentum trapped during the production phase for that event. 

The precise rules for the generation of the (") point are as follows. First, the C value 
for the point, ("*, is generated, according to a linear ramp distribution. This distribution 
extends between C = and C = 1; with value at ^ = and value 2 at C = 1. The ^ value 
for the point is then also generated. The distribution in this case is similar, except that 
now it extends between ^ = ^^(C*) and ^ = 1, ensuring that the point ends up inside the 
Yoshino-Rychkov bound. The details of how the program calculates ^6(C*) for the D and 
b appropriate to the event are given in Appendix B. 

The decision to implement a probability distribution favouring smaller mass losses 
than the Yoshino-Rychkov upper bound was made based on the results from the direct 
calculations given in Fig. 1. In this figure we have plotted the mean mass lost in an event 
with 6 = using the above probability distribution. The error bars represent the standard 
deviation in the b — mass loss. We observe a reasonably good agreement between the 
mean values obtained with our chosen probability distribution and the estimation method 
results, especially in the important D = 4 case where the estimation method results agree 
closely. Given that we favour smaller mass losses, it then seems sensible to ensure that the 
probability distribution also favours smaller angular momentum losses - hence the ramp 
distributions in ^ and (.^ 

^Recent results of simulations in four dimensions [59] indicate that ~ 25% of mass and ~ 65% of angular 
momentum are lost in collisions at the maximum impact parameter for black hole formation. However, 
the value obtained for the maximum impact parameter in this case is ~ 50% above the Yoshino-Rychkov 
lower bound, corresponding to an initial-state angular momentum that is more than double the maximum 
value possible for the black hole that is formed. Therefore an angular momentum loss greater than 50% is 
inevitable in this case. In Z3 > 5 dimensions there is no upper limit on the angular momentum of a black 
hole and such a large loss is not required. 
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Figure 2: Plots displaying the Yoshino-Rychkov bound (solid line) and some output C) points 
from our simulation of the production phase (dots), for selected D and b values (these are given 
above the plots in each case). Each plot contains 2000 sample output points, which have been 
generated with CVBIAS set to .TRUE.. 



One possible picture of the production phase is one in which the production phase 
radiation is 'flung out' radially in a frame co-rotating with the forming event horizon. 
In this scenario, the angular velocity of the event horizon does not change during the 
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production process (where we have to regard the forming black hole as a pseudo-Myers- 
Perry solution at all points during the production phase for this statement to have any 
meaning) . 

Based on this picture, we have included an option to bias the above probability dis- 
tribution such that C) points corresponding to smaller changes in the horizon angular 

velocity arc more likely. An additional condition is added to this bias — for the points to 
have their chance of being picked enhanced, they must also have an associated value of 
the oblatencss parameter a*, which is sufficiently close to that of the initial state. This 
is to remedy the problem that, for D > 5, there axe two curves with the same angular 
velocity as the initial state in the square 0<^<1,0<C<1, but only one is connected 
to the initial state. The bias may be turned on or off using the user switch CVBIAS, and 
the details of its implementation are discussed in Appendix C. 

In each of the plots in Fig. 2, there are 2000 C) points generated for the appropriate 
D and b using the new simulation with the bias applied. One can see in each case that 
there is an increased density of points around the 'constant horizon angular velocity' curve. 

The mass-energy lost during the production phase is distributed between radiation and 
the kinetic energy of the formed black hole. The production phase simulation must account 
for this. On the basis of several calculations [60-62], which indicate that gauge radiation is 
negligible compared to gravitational radiation in the production phase, we assume that all 
of the radiation is in the form of gravitons. Given that gravitons are missing energy, it is 
sufficient for the simulation to represent the entire radiation pattern using a 'net graviton' 
with a four-momentum equal to the sum of those of the emitted particles. 

The net graviton has an invariant mass /ig, which may potentially lie anywhere between 
and 1 — ^ (in units of the initial state mass) . An invariant mass of 1 — ^ corresponds to 
a completely symmetric emission of gravitons, whilst lower values correspond to steadily 
more antisymmetric emissions (which might result if a small number of gravitons is released, 
and by chance they are emitted in similar directions). In CHARYBDIS2, the invariant mass 
is randomly generated per event from a power distribution, P{iJ.g) oc Hg. The mean of this 
distribution is set equal to FMLOSTx(l-^) by the quantity FMLOST= {p+l)/{p+2) (default 
value 0.99, corresponding to p = 98). The simulation of the production phase emission is 
then a two body decay from initial state object into formed black hole plus net graviton, 
which is isotropic in the centre of mass frame of the initial state object. 

2.3 Adding the intrinsic spin of the colliding particles 

The results in the previous sections give a model for the angular momentum of the black 
hole after formation, which is based on using incoming particles with zero spin. Angular 
momentum conservation requires us to include the intrinsic spin of the incoming particles 
falling into the black hole. Since the results in the literature taking this effect into account 
are limited to special cases (see for example [51]) we assume a simple model where first we 
combine the spin states of the incoming particles into a state 

\si,hi)^®\s2,h2)^ = \s,Sz)^ . (2.2) 
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The collision axis is denoted by Sj, hi are the spin and helicity of the particles^ and s, Sz 
are the angular momentum quantum numbers of the combined state in the rest frame. 
Since we have unpolarised beams we give equal weight to each helicity combination. Then 
this angular momentum state is combined with the orbital contribution obtained from the 
model for angular momentum loss. We denote it by 

|L,L),, , (2.3) 

where L is the nearest integer to J. Note that z' is an axis in the plane perpendicular to the 
beam axis (z-axis) chosen with uniform probability. Finally, using a Wigner rotation [63] 
followed by a tensor product decomposition using Clebsch-Gordan coefficients, we obtain 



s 



s L+s 

= E^iU(O) E Cj,L,L,s,s,AJ,L + s,,,L,L,s,s,,)^, , {2A) 

s'^=—s J=\L—s\ 

(s) ■ ■ 

where dlj,^sz is a Wigner function and Cj^l^l^s,s^i is a Clebsch-Gordan coefficient for the 
tensor product decomposition of |L,L)^, (gi |s,s^/)^,. From (2.4) it is straightforward to 
determine the probabilities for all possible combinations of helicities and incoming partons. 

This model introduces a spread in the orientation of the initial black hole angular 
momentum axis around the plane perpendicular to the z-axis. Note that even though 
the model for angular momentum loss in the previous sections does not include such an 
effect, in a realistic situation we would not expect the angular momentum to be exactly 
perpendicular to the beam axis after the production phase. 

3. Decay of spinning black holes 

After the black hole settles down to the Myers-Perry solution, evaporation will start. The 
stationary background geometry for an (4 + n)-dimensional black hole with one angular 
momentum axis on the brane is [47] 

ds' =(l- ^) df + ^-^^^dtd<t. - ^dr'- 

_ ["^2 ^ ^2 ^ ^^S^) sin^ edcf" - cos^ Od^l , (3.1) 



where 



A = + ^ , S = + cos^ 9 , (3.2) 



t is a time coordinate, dJl^ is the metric on an n-sphere and {r, 6, (/>} are spatial spheroidal 
coordinates. This is clear if we transform to the coordinates {x, y, z} which define a spheroid 
(through their relation to r) according to 

4^ + 4 = 1 (3-3) 



*Note we are assuming the massless limit where hi = is. 
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at any space time point, and are expressed in terms of spheroidal coordinates in the asymp- 
totic r — > +00 region through 

X = Vr2 + a2 sin 9 cos ((/> + 0{^)) 

y = ^r'^ + a2 sin 9 sin {4> + 0{^)) . (3.4) 
z = rcos9 

The only two independent parameters in (3.1) are {/Lt, a} which are related to the physical 
parameters of the black hole. Those are respectively the extra dimensional generalisations 
of the ADM mass (M) and angular momentum (J)^ 



2 52+„(27r)-^M-+V, (3.5) 

n(n+l) „ 2 

J = 52+n(27r)-^J+^ M]5+2 a^i = — ^ Ma . (3.6) 

Furthermore, we can switch to a third pair of parameters, closely related to the geometrical 
properties of the black hole: {th, o*}- The first parameter is defined by the location of the 
horizon of the black hole at the largest positive root of A(rif) = 0. ru is directly related 
to the surface curvature of the horizon and thus (in some sense) has a frame independent 
meaning. The second is a* = ajru- This is easily interpreted as the oblateness of the 
spheroid at the horizon as seen from (3.3). 

We are working in an ADD scenario where, the Standard-Model fields are confined to 
a 4-dimensional brane. This type of setup has been proposed both for the case of fiat extra 
dimensions [1-3] (ADD scenario) and curved extra dimensions (for example the Randall- 
Sundrum models [4,5]). This constraint is necessary to avoid bounds from electroweak 
precision obscrvablcs. The upper bound on the extra dimensional width for such a brane 
is [64] R ~ (700 GeV)^^. This is many orders of magnitude below the scale needed in a 
large extra dimensions scenario (for n < 15) to explain the size of the Planck mass in the 
spirit of [2,3]. More elaborate setups have been proposed^ where different SM fields are 
placed on different branes imbedded in a thicker brane of size < 1 (TeV)~^ [65,66]. In 
our model we assume black holes of mass well above 1 TcV, which means that the typical 
Schwarzschild radius will be above 1 ~ 2(TeV)~^. Hence the minimum diameter should 
be 3 4 (TeV)^^, which is already well above the upper bound on the width of the thick 
brane. So compared to the size of the black hole, all branes should be well inside the black 
hole and for the purpose of evaporation, all fields effectively behave as being emitted from 
a single brane. Thus, for brane degrees of freedom, we use the 4-dimensional projected 
version of the metric (3.1), where the coordinates on the ri-sphere are fixed. Higher order 
corrections from splitting the branes are neglected. 

A related issue is that of brane tension. Here the results in the literature concerning 
transmission factors are limited to the case of a codimension-2 brane in six dimensions [67- 



^The expressions are obtained by looking at the metric when r — > +00 and comparing it with the metric 
for a weak localized massive perturbation in Minkowski space-time (see [47] for further details). 

®Such constructions aim to suppress the effects of dangerous operators which, for example, might allow 
fast proton decay. 
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71]. The main conclusion (see for example [67,69,70]) is that the brane projected metric 
(and consequently the Hawking spectra) remains unchanged up to a rescaling of the Planck 
mass, so only the spectra for bulk fields are qualitatively different. These observations, 
together with the fact that brane emission should be dominant, justify neglecting such 
effect for 6-dimensions (note that this case is anyway disfavoured from the bounds in [8,9]). 
For the phenomenologically favoured cases of codimension larger than 2, there arc virtually 
no detailed studies on brane tension effects. We may hope that a similar effect of rescaling 
of the Planck mass will occur for brane fields, in which case our model does not need to be 
adapted, but further work is required to understand if such an assumption holds. 

For our description of the evaporation to hold, we require the black hole to be sitting 
on a background that is effectively Minkowski spacetimc. This means that the typical size 
and curvature radius of the extra dimensions must be much larger than the horizon radius 
of the black hole. This holds for the ADD scenario, where the flat extra dimensions are 
large enough to solve the hierarchy problem, and for the Randall-Sundrum (RS) model 
with an infinite extra dimension [5], where the hierarchy problem is solved using a large 
curvature radius. 

We assume that the black hole remains stuck on the brane [72]. The possibility of 
ejection would come from graviton emission into the bulk [73,74]. Since the black hole is 
formed from SM particles which arc themselves confined, we assume that even if gravitons 
are emitted, the extra-dimensional recoiling momentum is absorbed by the brane, avoiding 
ejection. One could argue that whichever charges keep the black hole confined to the brane, 
they are lost at the start of the evaporation through Schwinger emission [16,17]. However 
it will be very unlikely that all the different gauge charges are simultaneously neutralized 
at any stage during the evaporation, if we assume that the black hole decays by emitting 
one quantum at a time. Furthermore, there are a lot more SM degrees of freedom than 
gravitational ones so even if the unlikely event of exact neutralisation occurs, it will still 
be unlikely that a graviton is emitted during the brief period of neutrality. Thus we would 
expect the number of events in which the black hole is ejected into the bulk to be at most 
a small fraction of the total. 

Even though exact neutralisation seems unlikely, it is well known, at least in four 
dimensions, that black holes tend to discharge very rapidly compared to the evaporation 
time. So we would expect the charge to stay low^, and charge effects on the probabilities 
of emission to be less important. In the generator we use a simplified model based on 
this observation, such that whenever a charged field is selected for emission, the electric 
charge of the state is selected as to reduce the total charge of the black hole (unless the BH 
is neutral, in which case equal probabilities for particles and anti-particles are used). To 
avoid complications in hadronization, baryon number conservation is assumed, and colours 
are assigned to ensure that colour singlet formation is possible. 

3.1 Energy distribution and greybody factors 

Once the evaporation starts, the black hole loses its mass and angular momentum through 
the emission of Hawking radiation. The radiation is thought to be predominantly in the 
'^For further references regarding discharge see chapter 10 of [80] and references therein. 
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form of SM fields on the brane, and gravitational radiation into the higher-dimensional 
bulk. It is believed [72] that "black holes radiate mainly on the brane", primarily due 
to the high multiplicity of the brane-confined SM fields. The conjecture is supported by 
studies of the ratio of power emission of scalar fields in brane and bulk channels [33,34]. 
Determining the energy balance between the brane and bulk channels for graviton emission 
remains a key open question. The emission of 'tensor' gravitational modes into the bulk 
was recently considered in [75,76]. A comprehensive analysis, including all gravitational 
modes ('tensor', 'vector' and 'scalar' [77-79]) for a simply rotating black hole, has yet to be 
conducted. In the following, we assume that emission on the brane is indeed the dominant 
decay channel. 

We assume further that the emission of Hawking radiation may be treated semi- 
classically. That is, that emission of radiation is a smooth, quasi-continuous process in 
which the black hole has time to reach equilibrium between emissions, and the energy of 
particles emitted is much less than the mass of the BH. These assumptions are valid as 
long as the mass of the black hole is much larger than the (higher-dimensional) Planck 
mass Md, in which case the typical time between emissions is large compared to rn (see 
for example (2.1)) and the Hawking temperature, which gives the typical energy scale of 
the emissions, is below Mr,. 

The Hawking temperature of a Myers-Perry black hole is 

(n+l) + (n-l)a^ 
= Ml + a^)r^ • ^'-'^ 

The particle flux, mass and angular momentum emitted by the BH per unit time and 
frequency in a single particle species of helicity h are 

d'{N,E,J} _ 1 ^ ' {l,u,m} 

dtdu - 2^ ^ ^ exp(^/r^)±l^^ ^^■^> 

j=\h\m=-j / -H; 

where A^, E, and J denote particle number, energy and angular momentum, respec- 
tively. Here, k = {h, j,m,uj} is shorthand for the numbers, and Co = lo — rn^l, where 
n = a*/[(l -I- a^)rif] is the horizon angular velocity. {h,j,m} are the spin weight, total 
angular momentum and azimuthal quantum numbers of the emission and u its energy. In 
the denominator, we select +1 for fermionic fields {\h\ = 1/2), and —1 for bosonic fields 
{\h\ = 0,1,2). The dimensionless quantity T^^^(a;,a*) is the transmission factor, also 
known as the greybody factor. 

In the context of Hawking radiation, a transmission or greybody factor is the proportion 
of the flux in a given mode that escapes from the horizon to infinity. No closed-form 
expression for these factors is known, although a number of useful approximations have 
been derived [30-32]. Simulation with CHARYBDIS2 requires accurate transmission factors 
across a wide parameter space; consequently a numerical approach was taken. 

To determine transmission factors one must solve the higher-dimensional Teukolsky 
equations [26,27,29]. After performing a separation of variables for a field on the brane. 
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the resulting radial equation is 




dRk 



dr 



) 



+ 



- ihKA' 
A 



(r) 



+ 4i/icc;r + /i(A"(r)-2),5;,,|;,|-Afe Rk = (3.9) 



where K = {r'^ + a'^)u) — am, and A = + — {r'jj + a^)r^^ /r"-~^ . Here K}^ = Ak — 
2mauj + a?' up', where Ak is the angular separation constant defined later on in (3.12). The 
transmission factors TT^^^ are found by considering modes which are purely ingoing at the 
outer horizon (for details see e.g. [26,27,29]). 

The separation constant is found by solving the angular equation (3.12) discussed in 
the next section. Unlike the radial equation, the angular equation does not depend on the 
dimensionality (£) = 4 + n) of the bulk spacetime. A range of methods for finding angular 
eigenvalues are detailed in [82]. For the purpose of calculation of transmission factors, we 
employed a spectral decomposition method (for example, see Appendix A in [83]) and a 
numerical shooting method [84]. We checked our results against known series expansions 
in aui [82]. 

The numerical methods employed to determine the transmission factors are described 
in detail in [26,27,29]. Transmission factors were computed numerically in the parameter 
range n = 1, 2, ... 6, ujth = 0.05, 0.10, ... 5.0 and a* = 0.0, 0.2, ... 5.0 for the angular modes 
j = \h\,\h\ + I, . . .\h\ +12 and m = —j . . . j. For each point we have computed the flux 
spectrum using (3.8). This quantity is used in CHARYBDIS2 as a probability distribution 
function for the quantum numbers of a particle with a given spin and to determine the 
relative probability of different spins (through integration of equation (3.8)). For conve- 
nience in the Monte Carlo, we have computed the following cumulative distributions from 
the transmission factors: 



where x is energy in units of r]^ , th = TffrH, x = x — mQ/rn and Q is an integer that 
counts modes. The modes are ordered with increasing j and within equal j modes they 
are ordered with increasing m. 

Cumulative functions are more convenient since they allow for high efficiency when 
selecting the quantum numbers. This is done by generating a random number in the range 
[0,C(oo)], followed by inversion of the corresponding cumulant. In CHARYBDIS2, when 
values of a* between those mentioned above are called, linear interpolation is used. When 
a* is larger than 5, we use the cumulative functions for a* = 5. We have checked that 
for most of the evaporation such large values are very unlikely. The exception is the final 
stage, when the black hole mass approaches the Planck mass. Here one of the remnant 
models takes over, as described in section 4. 

The effect of black hole rotation on the semi-classical Hawking emission spectrum is 
described in detail in the studies [23-29]. Here we briefly recall some qualitative features. 




(3.10) 
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(3.11) 
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Figure 3: Transmission factors and Planckian factors for a* = 0.8 and n = 6. The left plots show 
spin 1/2 and the right plots show spin 1. The top plots show the transmission factors and the 
bottom plots show the Planckian factors [exp((I'/T/f ) ±1] ^ , for a range of j and m modes. Lines 
with the same m have the same colour and line type. 

In the non-rotating (Schwarzschild) phase, the emission spectrum (3.8) is found from 
a sum over modes, and the modes are degenerate in m. In the 4D case, emission is 
dominated by the lowest angular mode {j = \h\), whereas in the higher-dimensional case, 
higher modes {j > \h\) are also significant. The Hawking temperature (3.7) increases 
monotonically with increasing n, and hence the total power (per particle species) increases 
steeply with increasing n. 

Rotation splits the azimuthal degeneracy; modes of different m are distinguished. Ra- 
diation is emitted preferentially into the co-rotating (m > 0) modes and causes the BH to 
lose angular momentum. The importance of co-rotating modes can be understood by exam- 
ining the interplay between the transmission factors (TFs) and the so-called Planckian 
(or thermal) factors {[exp{uj /Th) ± 1]^^) (PFs) in equations (3.8), as illustrated in Figs. 3 
and 4. At fixed j, we find the TFs for counter-rotating (m < 0) modes generally exceed the 
TFs for co-rotating modes, implying that (in some loose sense) the counter-rotating modes 
escape from the vicinity of the hole more easily that the co-rotating modes. However, the 
PFs are larger for co-rotating modes, implying that more Hawking radiation is generated 
in co-rotating modes. We find that this latter effect dominates over the former. In gen- 
eral, the effect of rotation is to substantially enhance emission, even though the Hawking 
temperature depends only weakly on a* for a fixed BH mass. 

For low dimensionalities (n = 0, 1, 2) and fast rotation (a* ~ 1) the emission spectrum 
is dominated by the maximally-corotating modes (m = +|j|). This can lead to an oscilla- 
tory or 'saw-tooth' emission spectrum, which is clearly visible in Fig. 4. The dominance of 
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Figure 4: Power spectra for a rotating six- dimensional black hole at a* = 1.4. The top plots show 
the transmission coefficients (red, solid) and the Planckian factors (blue, dotted) as a function of 
wr/i, for the m = j modes up to j = 19/2 for fermions (left) and j = 10 for vector bosons (right). 
The bottom plots display the power emission spectrum (thick red curve) containing contributions 
from all modes, together with the curves for the leading rii — j modes. The region of overlap between 
different modes is small, leading to sharply-peaked oscillations in the power emission spectrum. 

m = +\j\ modes will obviously lead to a rapid loss of angular momentum. If the number of 
bulk dimensions is large, n > 6, then this effect only occurs for very fast rotation a* ^ 1. 
In general for high n, all m modes contribute and the emission profile is much smoother. 

In the case of bosonic fields, black hole rotation induces another effect: the transmission 
factor T^^^ is negative for modes such that lolo < (see right-hand side plots of Figs. 3 
and 4) which can be checked from the Wronskian relations for the radial equation (3.9). 
That is, the transmitted part of 'in' wave modes with ujuj < falls into the rotating black 
hole carrying in negative energy and their reflected part returns to infinity with a gain in 
energy®. This classical phenomenon, which occurs for bosonic but not for fermionic fields, 
is known as superradiance [85-87]. Within the context of higher-dimensions, it has been 
shown [27] that spin-1 superradiance on the brane increases with the black hole intrinsic 
angular momentum a,, (as expected) and with the number n of extra dimensions; see 
also [25, 88] for studies of scalar superradiance on the brane. It has been suggested [89] 
that the presence of superradiance might lead to spin-2 radiation into the bulk being the 
dominant emission channel, thus disproving the claim that "black holes radiate mainly on 

*Note that the denominator in Eq. (3.8) is negative for superradiant modes in the bosonic case, thus 
cancelling out the negativity of the superradiant transmission factor T^^\ and so contributing positively 
to the fluxes. 
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the brane" made in [72]. However, it was shown in [34] (see also [90,91]) that superradiance 
for the scalar field emission in the bulk diminishes with the number n of extra dimensions, 
and that the percentage of the total (i.e., brane plus bulk) scalar power which is emitted 
into the bulk is always below 35% for n = 1, 2, . . . , 6. In the absence of a solution for the 
higher-dimensional spin-2 equations in the rotating case, these results for the scalar field 
seem to indicate that the main emission channel will be into the brane, rather than the 
bulk, in spite of superradiance. 

3.2 Angular distribution of Hawking radiation 

The angular distribution of brane fields emitted from a black hole rotating on the physical 
(3+l)-brane is controlled by the spheroidal wave functions Sk{c,x = cos6), which satisfy 
the differential equation 

d .2^rf^,„2.2 ("^ + ^^)' . - 5,(c,x) = 0, (3.12) 



J .{l-x^)—]+c^x^-2hcx-—_^ ^+Akic) + h 

ax \ ax J 1 — x'^ 

where c = auj and Ak{c) is the angular eigenvalue. In the generator, our method of solution 
for this equation is based on that of Leaver [92]: details may be found in Appendix D. 
Given the values of k and .4fc(c), the probability distribution function for an emission with 
momentum in the direction cos 9 is then given by the normalized square modulus^ of S^- 
This follows from the decomposition of spheroidal one-particle states into plane wave one- 
particle states, which is analogous to the decomposition of spherical waves into plane waves, 
for the usual case of scattering off a spherical potential. We work out the case of a massless 
Dirac field in Appendix E, in order to clarify the assignment of a spheroidal function of a 
certain spin weight to the correct helicity state in our convention. It turns out that in the 
convention of (3.12) the physical hehcity of the particle is actually —h. 

Some important properties of the angular distributions for our analysis follow from the 
observation of Fig. 5. We have verified that in general, for any mode, higher rotation tends 
to make the spheroidal functions more axial. This means that at low energies (where the 
modes are departing from being degenerate) , the angular distribution of Hawking radiation 
will tend to become more axial. However, for higher energies, the effect of rotation on the 
emission spectrum is to favour emission of modes with m = j in order to spin down the 
black hole. This will produce a more equatorial angular distribution, as we can see from 
Fig. 5 where the m = j mode is always more central (in x = cos 9) than the m = mode 
(for j / 0). So as we increase the rotation parameter we have a competition between the 
increase in the angular function's axial character and the increase in probability of emission 
of more equatorial modes (which are those with larger j). At low (high) energies the former 
(latter) wins. This observation is consistent with the energy dependence of the angular 
profiles shown in, e.g.. Fig. 16 of [29]. 

Low-energy vector bosons are more likely to be emitted close to the rotation axis, 
whereas high energy vector bosons are more likely to be emitted in the equatorial plane. 
A similar but far less pronounced effect exists for spin-half particles. 



'From now on, we will drop the explicit dependence on x and c whenever convenient. 
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Figure 5: The figure sliows various spheroidal wave functions Shj,m{x,auj) as a function of cos^. 
Plots for aoj = 0, 1, 2, 3, 4, 5, 6 are shown (only the first and last are indicated, since the curves are 
regularly ordered). The title of each plot indicates {h,j,m}. 



Particles with a single helicity, such as the neutrino, will be emitted asymmetrically 
by a rotating black hole [93]. For example, if the black hole's angular momentum vector 
is pointing north, then (anti-)neutrinos will be preferentially emitted in the (northern) 
southern hemisphere. W and Z boson decays will exhibit similar asymmetries in the two 
hemispheres. 
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3.3 Back-reaction and spin-down 

The difficult problem of studying the back-reaction is interesting, both from the theoretical 
and the phenomenological point of view, since it should start to influence the evaporation 
as the mass of the black hole is lowered. 

On the theory side, there is no well established framework to study the evolution of a 
Hawking-evaporating black hole over the full range of possible initial conditions. The usual 
approach [14, 94-96] is to write down mean value differential equations for the variation 
of the parameters (M and J) such as (3.8) and integrate them with appropriate initial 
conditions. However, this is only valid for a continuous process of emission where the 
variation of the parameters is very slow and the symmetry of the background spacetime is 
kept. Thus, it ignores the momentum recoil of the black hole and the change in orientation 
of the angular momentum axis between emissions, which are certainly negligible for an 
ultra-massive black hole, but will start to become important as we approach the Planck 
mass. Furthermore, since the Hawking spectra at fixed background parameters are used, 
it also neglects the effect of the backreaction on the metric by the emitted particle. This 
point has been explored in simplified cases of j = waves for fields of several spins using 
the method in [97] with some results regarding the modification of the thermal factors, but 
a full treatment is still lacking. 

In the program, we have included two possible models for the momentum recoil of the 
black hole set by the switch RECOIL, which takes the values 1 or 2. 

RECOIL = 1 interprets the selected energy as the energy of the particle in the rest 
frame of the initial black hole. The momentum orientation is computed in this frame with 
probability distribution given by the square modulus of the spheroidal function (3.12) and 
the momentum of the final black hole is worked out from conservation. The argument for 
this model comes from the observation that particles in the decay are highly relativistic. 
They propagate close to the speed of light, so the background they see is that of the initial 
black hole, since no signal of the back-reaction on the metric can propagate outwards faster 
than light. Thus the momentum of the emission is determined by the background metric 
in this picture. 

RECOIL = 2 takes the energy of the emission as being the loss in mass of the black 
hole. This corresponds to the usual prescription for computing the rate of mass loss. 
The orientation of the momentum in the rest frame of the initial black hole is computed 
as before with a probability distribution given by the square modulus of the spheroidal 
function (3.12) and the 4-momentum of the emission as well as that of the black hole are 
worked out. 

Note that for any of the previous options, full polarization information of the emission 
is kept, as it is generated with the correct angular distribution. This will potentially 
produce some observable angular asymmetries and correlations, which would not be present 
if angular distributions averaged over polarizations had been used. 

The other quantity we need to evolve is the angular momentum of the black hole. 
We have two options, controlled by the switch BHJVAR. The default BHJVAR = .TRUE, uses 
Clebsch-Gordan coefficients to combine the state of angular momentum = J of the 
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initial black hole (i.e. taking as quantisation axis the rotation axis of the black hole), 
with the emitted j, m state for the particle. The probability of a certain polar angle and 
magnitude for the angular momentum of the final black hole is given by the square modulus 
of the corresponding Clebsch-Gordan coefficient, and the azimuthal angle is chosen with 
uniform probability. If BH JVAR= . FALSE . , the orientation of the axis remains fixed even 
though the magnitude will change by subtraction of the m value of the emission. 

From previous versions of CHARYBDIS, we have kept the switch TIMVAR which allows 
one to fix the parameters of the black hole used in the spectrum (such as the Hawking 
temperature) throughout the evaporation. This option corresponds to a model where the 
evaporation is no longer slow enough for the black hole to re-equilibrate between emission, 
so in effect it represents a simultaneous emission of all the final state particles from the 
initial black hole without any intermediate states. 

In Fig. 6 we plot the evolution of the physical parameters M and J for BH events with 
fixed initial M, in the non-rotating case and the highly rotating case. In Fig. 7 we plot the 
horizon area and oblateness for the same cases as in Fig. 6. 

Events were generated for two possible initial masses, 10 TeV (reachable at the LHC) 
and 50 TeV. The latter serves as a check of the semi-classical limit. Wc focus on n = 2 
and n = 6. An important quantity necessary to produce these plots, is the time between 
emissions. Since our model for the evolution relics on the mean value equation (3.8), before 
each emission, an average time can be computed (i.e. St for 5N = 1): 

-1 



SN = —St ^ St = 
dt 



dN 



where a sum over all species is assumed (see Sect. 4 for further details). 

Each plot contains 10^ trajectories (one per event generated), each contributing with 
weight 1 to the density plot. The darker areas correspond to higher probability and in all 
the plots we can discern a tendency line which is sharper for the 50 TeV case and more 
diffuse for 10 TeV, reflecting the magnitude of the statistical fluctuations. 

The left columns of Figs. 6 and 7 show respectively the evolution of the mass parameter 
and horizon area for non-rotating black holes. The centre and right columns show the 
evolution of the mass and angular momentum, or the horizon area and oblateness, for the 
highly rotating case (a* ~ 3). The main features are as follows: 

• Non-rotating case: Both M and A decrease approximately linearly with time except 
for the last ~ 10 — 20 % when they drop faster. This is directed related to the 
behaviour of the temperature which increases slowly (approximately linearly) for 
most of the evaporation and rises sharply near the end. The rates tend to be faster 
for higher n which is in agreement with the increase in Hawking temperature with n. 



Highly-rotating case: Here the statistical fluctuations tend to smear out the plots 

for the case of lowest mass. However, the same tendency can be seen as for the 
M = 50 TeV black holes; the latter display better a true semi-classical behaviour. 
There is in general an initial period of roughly 10 — 15 % of the total time when M 
drops faster to about 60 — 70%. At the same time the angular momentum also drops 
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Figure 6: Probability maps for pliysical parameters, constructed from 10^ trajectories for different 
BH events with fixed initial conditions M, J (for each horizontal line). Each trajectory contributes 
with weight 1 to the bins it crosses on the {P,t/ttotai} plane where P is the relevant parameter. 
Note that the time is normalised to the total time for evaporation ttotai- The horizontal lines for 
the plots on the right are due to the discretisation of J in semi-integers. 



sharply to 20 %. This corresponds to the usual spin-down phase [95]. Note that 
the fluctuations are quite large for the low-mass n = 2 plots. As for the geometrical 
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Figure 7: Probability maps for the geometrical parameters A and a* which characterize respec- 
tively the size and oblateness of the BH. Note that for the a* plot in the last t/ttotai = 1 line there 
is very often a jump to very large a,. In these plots we have put all such points in the bin on the 
upper right corner to avoid squashing the interesting region. Each horizontal line has the same 
initial M, J as the corresponding one in Fig. 6. 



parameters, they follow a similar tendency if we make the correspondences M ^ A 
and J ^ a^. Again for the low- mass plots the statistical fluctuations smear out 
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the sharper initial drop in area, in particular for n = 2 in which it can occasionally 
increase substantially. The a* plots show how the black hole tends to become more 
spherical in this spin-down phase. The remainder of the evolution resembles the non- 
rotating case and can be identified with a Schwarzschild phase. Note however that 
by the end of the evaporation (when M approaches the Planck mass) this description 
breaks down and rises again, since even only one unit of angular momentum has 
a very large effect on this quantity at the Planck scale. This means that a* ceases 
to have a well defined geometrical meaning, as we reach the Planck phase. Similarly 
to the non-rotating case, we have checked that the temperature increases slowly and 
approximately linearly for most of the evaporation except for a sharp rise as the 
Planck mass is approached. 

These observations agree with the usual results in four dimensions'^ and the results 
of [28] in D dimensions. 

4. Termination of black hole decay 

Our model for black hole decay relies heavily on the assumptions that we are in the semi- 
classical regime and the evaporation is slow (i.e. there is enough time for rc-cqiiilibration 
between emissions) [6,98]. However, as the evaporation evolves, we will reach a point where 
neither of these assumptions will be true. 

In the generator we introduce some options for the final remnant decay based on 
different physical assumptions. First of all, we need a criterion to decide whether or not 
the remnant stage has been reached. The various options in the program are connected 
to a departure from semi-classicality. This occurs when the expectation value (N) for the 
number of emissions becomes small, which is a sign of the low number of degrees of freedom 
associated with the black hole. Together with the drop in (N), the Hawking temperature 
will rise sharply. This is all related to the approach of the black hole mass to the Planck 
mass. The options are: 

• NBODYAVERAGE= . TRUE . : An estimate for the multiplicity of the final state is com- 
puted at each step during the evaporation, according to the Hawking spectrum: 



The sums are over all particle species with appropriate degeneracies gi. The inte- 
grated flux and power are computed using (3.8). A natural criterion for stopping 
the evaporation is when this estimate drops below some number close to 1. In the 
generator we use {N) < NBODY — 1 where NBODY gives the average multiplicity of the 
remnant decay final state (see Sect. 4.2 for further comments). Varying the param- 
eter NBODY will give a measure of uncertainties in the remnant model. In addition, 

see for example [95] or chapter 10.5.3 of [80] 




(4.1) 
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if wc choose a remnant model that decays, (4.1) gives an estimate of the final state 
multiplicity for such a decay. When NBODY AVERAGE = .FALSE., one of the options 
below, inherited from earlier versions of CHARYBDIS, is used. 

• KINCUT= . TRUE . : Terminate evaporation if an emission is selected which is not kine- 
matically allowed. This is closely related to the rapid increase in temperature as we 
approach the Planck mass and consequently the generation of kinematically disal- 
lowed energies for the emission. Otherwise if KINCUT = .FALSE, the kinematically 
disallowed emissions are rejected and the evaporation terminates when the mass of 
the black hole drops below the Planck mass^^. 



4.1 Fixed- multiplicity model 

The default remnant decay option is a fixed-multiplicity model similar to that in earlier 
versions of CHARYBDIS. At the end of the BH evaporation, the remaining object is decayed 
isotropically in its rest frame, into a fixed number NBODY of primary particles, where the 
parameter NBODY is an integer between 2 and 5. The decay products are chosen with 
relative probabilities appropriate to the final characteristics of the black hole (i.e. weighted 
according to the integrated Hawking fluxes for each spin). 

The selection of the outgoing momenta of the decay products may be chosen either 
using pure phase space (NBODYPHASE= . TRUE . ) or by using the following probability density 
function in the rest frame of the black hole (NBODYPHASE=. FALSE.): 

dP oc (5(^) (^pi - Pbh^ n Pi (-^^' d^^P^ ' (^-2) 

which amounts to the usual phase space momentum conservation with an extra weight 
function for each particle 

p, (E,,a) = \Sk{cos9i)\' , (4.3) 

exp[E/TH) ± 1 

where k = {j,m} are chosen according to the cumulants (3.11) combined with angular 
momentum conservation. Here Ei,Qi are the energy and momentum orientation of the 
emission in the rest frame of the remnant. This choice treats the final state particles on 
an equal footing, keeping a gravitational character for the decay (since it uses Hawking 
spectra) , as well as some correlations with the the axis of rotation through the spheroidal 
function factor. Furthermore, at this stage, slow evaporation should no longer be valid, 
so it makes sense to perform a simultaneous decay at fixed black hole parameters. This 
remnant option can be used with any of the criteria for termination. 



^^The Planck mass used in CHARYBDIS2 to decide on the termination is always the internal one, INTMPL, 
which is obtained by converting the Planck mass input by the user (in a given convention), to the Giddings- 
Thomas convention - see Appendix A. 
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4.2 Variable- multiplicity model 

In addition to a fixed multiplicity final state, an option has been introduced to select the 
multiplicity of the final state on an event-by-event basis. Wc follow an idea in [99], which 
has been used for example in the case of 2 ^ 2 sub-processes in [44]. Here we implement 
a more general model for arbitrary multiplicity, which is invoked by setting the parameter 
NBODYVAR=.TRUE.. 

As argued previously, when the remnant stage is reached, the black hole should no 

longer have time to re-equilibrate between emissions. Under this assumption, the prob- 
ability distributions should become time independent. It is relatively straightforward to 
prove that under these conditions for a time interval 6t, the multiplicity follows a Poisson 
distribution [99]: 

Pstin) = e--''^^ , (4.4) 
n! 

with a some constant. From the Hawking flux, we have computed an estimate for the aver- 
age number of particles emitted during 6t (i.e. the time interval until all mass disappears), 
so a is determined from this condition. The final result is 

P.t(n) = e-Wi^ ' (^-^^ 

where {N) is the estimate in (4.1). This expression gives us an estimate for the probability 
of emission of n particles from the remnant, so we choose to interpret n+1 as the multiplicity 

of the final system. In the generator we have removed the n = case (i.e. multiplicity 1 
final state) since the probability of the remnant to have all the correct quantum numbers 
and mass of a standard model particle will be vanishingly small. 

After the multiplicity is chosen, cither the pure phase space decay or the model de- 
scribed in the previous section is used, according to the value of NBODYPHASE. 

4.3 Boiling model 

The boiling remnant model, activated by setting RMBOIL= .TRUE. , is loosely motivated by 
the expectation that at the Planck scale the system becomes like a string ball [41,100], which 
has a limiting temperature due to the exponential degeneracy of the string spectrum [45]. 
In this model, evaporation of the BH proceeds until the Hawking temperature for the next 
emission would exceed a maximum value set by the parameter THWMAX. Prom that point on, 
the temperature is reset to THWMAX and the oblateness is frozen at the current value. The 
remaining object evaporates like a BH with those characteristics, until its mass falls below 
a value set by the parameter RMMINM. It then decays into a fixed number NEDDY of primary 
particles, as in the fixed-multiplicity model, or a variable number if the variable-multiplicity 
model is on. 

4.4 Stable remnant model 

A number of authors have proposed that the endpoint of black hole evaporation could 
be a stable remnant [42,43,46]. This option is activated by setting RMSTAB= . TRUE . . In 
order for the cluster hadronisation model of HERWIG to hadronise the rest of the final state 
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successfully, the stable remnant must be a colourless object essentially equivalent to a 
quark-antiquark bound state. Therefore it is required to have baryon number Br = and 
charge Qr = or ±1. 

The stable remnant appears in the event record as RemnantO, Remnant+ or Remnant-, 
with PDG identity code 50, 51 or -51, respectively, according to its charge. This object 
will behave as a heavy fundamental particle with conventional interactions in the detector. 

If a remnant with Br 7^ or \Qr\ > 1 is generated, the whole BH evaporation is 
repeated until Br = and \Qr\ < 1. This can make the stable remnant option much 
slower than the other options, depending on the length of the black hole decay chain. 

4.5 Straight-to-remnant option 

Recently, there has been discussion of the possibility that the formation of a semi-classical 
black hole may lie beyond current experimental reach, with low-multiplicity gravitational 
scattering more likely at the TeV scale [44]. To simiilatc this scenario, CHARYBDIS2 provides 
the option of bypassing the evaporation phase by setting the switch SKIP2REMNANT= . TRUE . 
and skipping directly to one of the remnant models presented in the previous sections. This 
permits the study of a wide range of qualitatively different possibilities, from simple 2^2 
isotropic scattering (fixed multiplicity) to more complicated variable-multiplicity 2 ^ N 
sub-processes. 

The 2 — > iV model is particularly flexible, allowing cither a phase-space distribution 
or one using the Hawking energy and angular spectra (see Sect. 4.2). Apart from this, 
all particle species are treated on an equal footing consistent with conservation laws. Al- 
ternatively the quantum-gravity motivated boiling model can be used. Further work will 
be presented in future publications exploring the phenomenological consequences of these 
scenarios. 

5. Results 

In this section we present results from CHARYBDIS2 simulations of black hole production at 
the LHC. A range of CHARYBDIS2 samples were produced using HERWIG 6.510 [101, 102] to 
do the parton showering, hadronisation and standard model particle decays. The results 
of which were then passed through a generic LHC detector simulation, AcerDET 1 . [103]. 
CHARYBDIS2 parameter defaults are shown in Table 1. In all following discussion, the 
number of extra spatial dimensions is n = TOTDIM — 4. Samples were generated with 
a 1 TeV Planck mass (in the PDG convention, i.e. MSSDEF = 3) so as to investigate the 
phenomenologically preferred region accessible at the LHC. Black holes were generated with 
a lower mass threshold of 5 TeV such that the semi-classical approximations for production 
are valid. 

Our settings for AcerDET 1 . are as follows: we select electrons and muons with 
Pt > 15 GcV and \rj\ < 2.5. They are considered isolated if they lie at a distance 

AR = \J (A?7)2 -|- (A(/>)^ > 0.4 from other leptons or jets and if less than 10 CeV of energy 
was deposited in a cone of Ai? = 0.2 around the central cluster. The same prescription 
is followed for photons. Jets are reconstructed from clusters using a cone algorithm of 
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Table 1: Default CHARYBDIS2 generator parameters (new parameters in the second set). 



Name 


Description 


Default 


MINMSS 


Minimum parton-parton invariant mass 


5 TeV 


MAXMSS 


Maximum parton-parton invariant mass 


14 TeV 


MPLMCK 


Planck scale 


1 TeV 


GTSCA 


Use Giddings-Thomas scale for PDFs 


.FALSE. 


MSSDEF 


Convention for Planck scale 


3 


TDTDIM 


Total number of dimensions 


6 


NBODY 


Number of particles in remnant decay 


2 


TIMVAR 


Allow Th to evolve with BH parameters 


. TRUE . 


MSSDEC 


Allowed decay products (3=all SM) 


3 


GRYBDY 


Include grey-body factors 


.TRUE. 


KINCUT 


Use a kinematic cut-off on the decay 


.FALSE. 


THWMAX 


Maximum Hawking temperature 


1 TeV 


BHSPIN 


Simulate rotating black holes 


. TRUE . 


BHJVAR 


Allow black hole spin axis to vary 


. TRUE . 


BHANIS 


Non-uniform angular functions for the evaporation 


.TRUE. 


RECOIL 


Recoil model for evaporation 


2 


MJLOST 


Simulation of M, J lost in production/balding 


.TRUE. 


CVBIAS 


'Constant angular velocity' bias 


.FALSE. 


FMLDST 


Isotropy of gravitational radiation lost 


0.99 


YRCSC 


Use Yoschino-Rychov cross-section enhancement 


. TRUE . 


RMSTAB 


Stable remnant model 


. FALSE . 


NBODYAVERAGE 


Use flux criterion for remnant - see Eq. (4.1) 


.TRUE. 


NBODYVAR 


Variable-multiplicity remnant model 


.FALSE. 


NBODYPHASE 


Use phase space for remnants 


.FALSE. 


SKIP2REMNANT 


Bypass evaporation phase 


.FALSE. 


RMBOIL 


Use boiling remnant model 


. FALSE . 


RMMINM 


Minimum mass for boiling model 


100 GeV 



AR = 0.4, with a lower Pt cut of 20 GeV. Lepton momentum resolutions were parame- 
terised from ATLAS full simulation results published in [104].^^ Where reference is made 
to reconstructed multiplicities or spectra, the reconstructed objects are either electrons, 
muons, photons or jets from AcerDET. 

5.1 Black hole mass and angular momentum 

The cross-section for black hole production is a strong function of the Planck mass. Though 
not affecting the total black hole cross-section, simulating the mass and spin lost during 
black hole formation does have a large effect on the cross-section for a particular mass 

Electrons are smeared according to a pseudorapidity dependent parameterisation; for muons, we take 
the resolutions from \ri\ < 1.1. 
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Figure 8: Differential cross-sections for black hole production with n = 4 extra dimensions. They 
are shown for different settings of the Planck mass and MJLDST (simulation of mass and angular 
momentum lost in production/balding). 

range. The differential cross-section will be reduced, for the same input state will produce 
a black hole of lesser mass, as is illustrated in Fig. 8. 

Most collisions with sufficient energy to create a black hole are between two (valence) 
quarks, however a minority occur in collisions between a quark and a gluon. CHARYBDIS2 
adds the spins of the colliding partons when forming a black hole; the initial angular 
momentum is either integer or half-integer accordingly. An integer loss of orbital angular 
momentum in the formation process is simulated by the Yoshino-Rychkov model described 
in Sect. 2. 

At high n there is a large increase in the production of high spin states, as seen in 
Fig. 9. The average spin of the produced black hole rises from 5.0 units for n = 2, to 8.1 
for n = 4 and 10.6 for n = 6. The generator switch MJLOST toggles a model of the loss of 
black hole mass and angular momentum in graviton emission during black hole production 
and balding, as described in Sect. 2.2. Setting this to .TRUE, decreases the spin slightly by 
an average of 30% for n = 2, 4, 6, whilst the mass drops by 18% (n = 2) to 30% (n = 6), 
as shown in Fig. 9. 

The variation in the number of Hawking emissions is caused primarily by the black 
hole mass, though the black hole spin and temperature play a role - a more highly rotating, 
or higher temperature black hole will emit more energetically. Consequently, the decrease 
in the number of Hawking emissions follows the drop in mass and is greatest for higher 
numbers of extra dimensions, with an average of two fewer emissions (or 30%) manifest for 
n = 6. 

The simulation of losses from production and balding has several major effects upon 
the produced particle spectrum. The reduction in the black hole mass and number of 
Hawking emissions leads to a decrease in the number of particles observed experimentally 
and to a reduced differentiation between samples with different numbers of dimensions. 
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Figure 9: Effects of mass and angular momentum loss in the formation/balding phase. The left 
(right) upper plot shows the angular momentum distribution before (after) this phase, whilst the 
lower row details the resulting black hole masses and the number of Hawking emissions. 



The decrease in initial black hole mass also leads to a softening of the emitted particle 
spectrum, with the high energy and transverse momentum tail of the distribution being 
reduced. Fig. 10 shows this for a fixed 2-body remnant decay using the option in the second 
row of Table 2. 



Table 2: Parameters used for remnant comparison. All samples have n = 2 and MJLOST=. FALSE. 



Legend 


Remnant Criterion 


Fixed /Variable 


Remnant No. /Mean 


Kincut on 


M < INTMPL (KINCUT= . TRUE . ) 


Fixed 


2 


Kincut off 


M < INTMPL (KINCUT=. FALSE.) 


Fixed 


2 


Nbody2 


Flux (NBODYAVERAGE= . TRUE . ) 


Fixed 


2 


Nbody3 


Flux (NBODYAVERAGE= . TRUE . ) 


Fixed 


3 


Nbody4 


Flux (NBODYAVERAGE= . TRUE . ) 


Fixed 


4 


Nvar2 


Flux (NBODYAVERAGE= . TRUE . ) 


Variable 


2 


NvarS 


Flux (NBODYAVERAGE= . TRUE . ) 


Variable 


3 


Nvar4 


Flux (NBODYAVERAGE= . TRUE . ) 


Variable 


4 


Boiling 


RMMINM < M < INTMPL 


Variable 


2 
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Figure 10: Effect of simulating tlie mass and angular momentum lost in blaek hole production 
on particle multiplicity distributions and spectra at generator level (left) and after AcerDET 
detector simulation (right), for a fixed 2-body remnant decay using the option in the second row of 
Table 2. 

5.2 Rotation effects 

The inclusion of black hole angular momentum has several large effects upon the emitted 
particles and their spectrum. The probability of a highly energetic emission is enhanced for 
partial waves with high values of the azimuthal quantum number m. As discussed at the 
end of Sect. 3.1, this results from the interplay between Planckian factors and transmission 
coefficients. The former, which turn out to be the dominant effect, are enhanced strongly 
for large positive values of m due to the —mil. term in the exponential of equation (3.8), 
reducing the Planckian suppression. Consequently, the particle energy and transverse 
momentum (Pt) distributions for emissions from a rotating black hole are harder. The 
number of primary emissions is correspondingly reduced. Fig. 11 shows the emitted particle 
multiplicity and Pt spectra for different numbers of extra dimensions. The effects of black 
hole rotation are largest for fewest number of extra dimensions, for which the spin term 
(in the Planckian factors) has greater magnitude. This more than compensates for their 
slightly lower Hawking temperature. 

The effect of black hole rotation on the pseudorapidity (rj) distribution is more subtle. 
Assuming no strong spin recoil during the balding phase, the initial black hole formed will 
have a spin axis perpendicular to the beam direction. 
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Figure 11: Particle multiplicity distributions and Pt spectra at generator level (left) and after 
AcerDET detector simulation (right) for non-rotating and spinning black hole samples, with rt — 2, 
4 and 6 extra dimensions and MJLOST = .FALSE.. 




Figure 12: Normalised particle rj distributions at generator level (left) and after AcerDET detector 
simulation (right) from black hole samples with n extra dimensions and MJLOST = .FALSE.. 



Since emission in the equatorial plane is favoured, particularly for scalars and fermions, 
one would expect the component along the beam direction, and hence at high rj, to be 
enhanced, at least for initial emissions. This effect is seen experimentally in (Fig. 12), but 
is slight. 

Similar trends can be seen in event variables such as missing transverse energy (MET) 
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Figure 13: Missing transverse energy and scalar Pt sum for rotating and non-rotating black hole 
samples after AcerDET fast detector simulation. Samples used the NBODYAVERAGE criterion for the 
remnant phase and include a simulation of the mass and angular momentum lost during production 
and balding. 



and The reduced particle multiplicity increases the probability of minimal or no 

MET, where no neutrinos are present in the event (neither directly emitted by the black 
hole, nor in weak decays of other primary emissions). The greater energy of the Hawking 
emissions increase the very high MET tail: a neutrino emitted by a rotating black hole 
is likely to have higher energy and momentum. The result is a flatter, longer tail for the 
rotating case, extending further beyond 1 TeV, as shown in Fig. 13. 

When compared to the number of primary emissions from the evaporation, a greater 
number of detector objects (leptons, photons, hadronic jets) are observed following fast 
detector simulation. Neutrinos emitted by the black hole will not be seen experimentally, 
whereas a single heavy quark or vector boson will result in the detection of multiple particles 
or jets of hadrons. Equally, the transverse momentum spectrum observed experimentally 
will be slightly softer in general than that of the primary particles emitted by the black 
hole, due to secondary emissions and radiation. 

5.3 Particle production probabilities 

Black hole rotation has a large effect on the particle production probabilities (Fig. 14). 
The most dramatic of these is the enhanced emission coefficient for vector particles. This 
is due to the larger fluxes and agrees with the greater differential power fluxes per degree 
of freedom shown in Fig. 4 (see vertical axis). 

The greater proportion of vector emissions would provide strong evidence of rotating, 
rather than Schwarzschild, black holes. However such measurements are difficult to make 
in practice - at the LHC it will not be possible to distinguish gluon jets from quark ones. 
Though highly boosted vector bosons provide experimental challenges, Z bosons can often 
be studied via their leptonic decay modes. Perhaps the most accessible other means to 
investigate black hole rotation might be the study of the photon multiplicity or its ratio to 
other particles, TeV-energy photons being one manifestation of black holes reproduced by 
neither other new physics scenarios nor SM backgrounds. Another experimental difficulty 
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Figure 14: Particle emission probabilities for rotating and non-rotating black holes with n — A 
(top) and for rotating black holes with different n (bottom) on linear and logarithmic scales. The 
lower left figure shows the fractional enhancement of specific particle emission from rotating black 
holes, as compared to non-rotating ones. The x-axis shows the PDG identification code for the SM 
particles, as defined in [f05]. 



for the detection and the isolation of the black-hole signal is that rotation decreases the 
probability of producing a lepton - often useful in reducing jet-like SM backgrounds to 
black hole events [106]. 

The emission probabilities for each particle species are largely independent of the 
number of extra dimensions, which primarily affects the emission energy and multiplicity, 
so that a reproduction of the distribution of particle species would be powerful evidence of 
black holes. 

The particle-antiparticle imbalance in Fig. 14 is chiefly caused by the (usually positively 
charged) input state. According to the mechanism described in Sect. 3, up-type quarks 
and down-type antiquarks are favoured, so as to meet the constraints of charge balance. 
Similarly, the net positive baryon number of the input state and the need to conserve 
baryon number for hadronisation leads to a preference for quarks over antiquarks. The 
apparent increase in this with rotation is a reflection of the reduced particle multiplicity: 
with fewer particles amongst which to share the charge imbalance, the effect is magnified. 
This is potentially a source of uncertainty since, unlike charge, black holes do not have to 
conserve baryon or lepton quantum numbers. At present we are constrained to conserve 
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Figure 15: The evolution of black hole parameters during the Hawking evaporation phase. The 
lower two rows show the mean value tendency curves of the relevant parameter as a function of 
fractional time t/tTotai through the evaporation. Error bars indicate the standard deviation of the 
distribution. The bottom right plot indicates the typical time interval 5t between emissions; its 
large standard deviation is due to a long tail for the distributions at each t/tTotai- 



these by the needs of hadronisation generators. 
5.4 Black hole evolution 

As the Hawking emission proceeds, the black hole evolves, becoming lighter, hotter and 
losing angular momentum as detailed in Fig. 15. Without the simulation of losses in pro- 
duction/balding, the black hole spins down more quickly than it loses mass; for large black 
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Figure 16: 2D contour plots showing how the black hole mass, angular momentum evolve with 
successive Hawking emissions, and their correlations for samples with n — A extra dimensions, 
using the NBDDYAVERAGE criterion for the remnant phase and including a simulation of the mass 
and angular momentum lost during production and balding. 

hole angular momentum, energetic emissions at high m are highly favoured. Turning on 
the simulation suppresses initial states with high black hole angular momentum. Conse- 
quently this effect is reduced in magnitude, though the majority of the black hole angular 
momentum is still lost before its mass. 

Initially the high angular momentum of the black hole leads to a high emission flux so 
the typical time interval between emissions is reduced relative to a non-rotating BH. 

The higher spin term in the Planckian factor causes there to be fewer, more energetic 
emissions for few extra dimensions. Half the mass is lost in the first 3 emissions for n = 2, 
compared with 4 (n = 4) and 5 for n = 6. The distribution does have a substantial tail 
however, with 1% of black hole events producing more than 11 primary emissions. 

As is shown in the contour plots of Fig. 16, black holes with high initial angular 
momentum tend to lose much of it during their first few emissions, whereafter further 
emissions decrease the black hole mass more smoothly, whilst its angular momentum stays 
relatively low, but non-zero. Thus the spin-down phase persists throughout the black hole 
decay - only a small proportion of black holes settling into a Schwarzschild, non-rotating 
state. This is in direct agreement with the theoretical plots in Figs. 6 and 7. 

As the black hole becomes lighter, its temperature rises, as does its oblateness (a*) 
and the time interval between emissions drops. These effects are gradual except when the 
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Figure 17: Primary particle multiplicity and \Pt \ distributions for black hole samples with n — 2, 
using a wide range of remnant options, as defined in Table 2. 



black hole mass becomes very low, at the end of the Hawking radiation phase. At this 
point we have reached the remnant phase. 

The striated, lined structure in the angular momentum plots is due to the initial state: 
the probability of a quark-gluon interaction, and correspondingly a half-integer angular 
momentum state is much lower than the integer state above and below it. 

5.5 Remnant options 

CHARYBDIS2 includes several models for the remnant phase. Both fixed and variable body 
decays have parameter switches to enable the systematics to be studied, as detailed in 
Sect. 4. 

The fixed multiplicity model, present in CHARYBDIS and optional in CHARYBDIS2, is 
linked to the choice of the variable KINCUT. If KINCUT=. FALSE. , proposed decays that 
are kinematically disallowed are ignored; if KINCUT= . TRUE . , their proposal terminates the 
evaporation phase. The former choice will give a greater number of less energetic particles, 
as evidenced by Fig. 17 which constrasts a range of remnant models defined in Table 2. 

CHARYBDIS2 uses the NBODYAVERAGE remnant criterion as a default, where the fluxes 
are used to calculate the expected number of further emissions. This provides a physically 
motivated model. Using this criterion with either a fixed 2-body ( "nbody2" ) or a variable 
multiplicity remnant model ("nvar2") gives a distribution lying between the upper and 
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Figure 18: Particle type and Pt plots for 2-body remnant decays using the old model "Kincut off" 
(top) and the new model "nbody2" (bottom) as defined in Table 2. Distributions are normalised 
per event. 



lower values obtained using the older model switches (upper plots of Fig. 17), indicating 
good control over the uncertainties mentioned in Sect. 4. The string-motivated boiling 
model gives a slightly higher multiplicity, since successive emissions are produced until the 
remnant mass drops below the remnant minimum mass, resulting in a greater number of 
softer particles produced in the remnant phase. 

The new NBODYAVERAGE model is also more robust with respect to changes in the 
number of particles produced in the remnant phase. This is because the flux calculation 
allows the spin-down phase to be terminated whenever the expected number of further 
emissions is fewer than that selected for the remnant phase. This is illustrated in the lower 
plots of Fig. 17, where changing the number of particles produced in the remnant phase 
results in similar multiplicities and spectra; events with 4-body remnant decays do not 
always have two more particles than their 2-body analogues. 

Another advantage of the NBODYAVERAGE method is that by using the integrated power 
and flux, the spin-down phase is terminated at a point that allows a smoother transition to 
the remnant phase, as shown by their spectra in Fig. 18, where the NBODYAVERAGE (lower) 
method gives a more concordant distribution of particle transverse momenta. Performing 
a remnant decay only when the mass drops below the Planck mass gives a much softer 
momentum spectrum, in contrast to the high energies favoured by light rotating black 



-36- 




Reco Mass - True Mass [GeV] Reco Mass - True Mass [GeV] 

Figure 19: Sample black hole mass resolutions after AcerDET detector simulation with n — 2 and 
no balding simulation for all events (left) and after a cut of MET< 100 GeV (right). The fits are 
indicative of the resolution in the peak and do not model the non-Gaussian tails which remain. 



holes. The option to start the remnant decay based on the drop of (iV) provides a smoother 
transition, since the final decay particles will have a harder spectrum, more similar to the 
Hawking phase. Emissions in the remnant phase are predominantly coloured, with positive 
baryon number favoured, so as to meet the constraints of baryon number conservation. 

5.6 Mass reconstruction 

In principle, it is possible to reconstruct the black hole mass by combining the 4-momenta 
of all particles observed in the event and missing transverse energy. Mass resolutions of 
200-300 GeV can be achieved for some samples as shown in Fig. 19, but there is significant 
variation with different samples and black hole parameters. Events with large amounts of 
MET (particularly from multiple sources) tend to be more poorly reconstructed. Invoking 
a 100 GeV cut on MET results in better reconstruction at the cost of some signal efficiency. 
Such a cut may not be entirely conservative however, for there may be additional sources 
of MET neglected in our simulation, such as that from Hawking emission of gravitons. 

5.7 Planck mass conventions 

Although a one-to-one mapping between differing conventions for the Planck scale is 
straightforward to make, naively changing the convention without compensating for the 
Planck mass value leads to large apparent differences in the output distribution, due to their 
different n dependence. This is particularly the case for high numbers of extra dimensions 
for which the definitions vary most widely, as shown in Fig. 20 where the Dimopoulos- 
Landsberg and PDG Planck mass definitions are contrasted. 

A 1 TeV Planck mass in the Dimopoulos-Landsberg convention equates to a larger 
value in the PDG convention, especially for large n (the increase ranges from a factor of 
1.1 at n = 1 to 2.7 at n = 6). Consequently they have a higher Hawking temperature and 
far fewer particles are emitted, but with higher energies and transverse momenta. 
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Figure 20: Particle Multiplicity distributions and Pt spectra at generator level (left) and after 
AcerDET detector simulation (right) for a 1 TeV Planck mass in PDG and Dimipoulos-Landsberg 
"MDL" conventions. 



6. Conclusions 



We have presented in detail the physics content of the new black hole event generator 
CHARYBDIS2, together with some results illustrating important features of the simulation 
of the different phases of black hole production and decay. The main new features compared 
to most earlier generators, including CHARYBDIS, are: detailed modelling of the cross section 
and the loss of energy and angular momentum during formation of the black hole (the so- 
called balding phase), based on the best available theoretical information; full treatment of 
angular momentum during the evaporation phase, including spin of the incoming partons, 
rotation of the black hole, and anisotropy and polarisation of all Standard-Model fields 
emitted on the brane; and finally a variety of options for the Planck-scale termination 
phase, ranging from a stable remnant to a variable-multiplicity model connecting smoothly 
with the evaporation phase. 

Our main finding is that angular momentum has strong effects on the properties of the 
final state particles in black-hole events. Even after allowing for a substantial loss of angular 
momentum in the balding phase, the isotropic evaporation of a spinless Schwarzschild- 
Tangherlini black hole is not a good approximation, nor is the notion of a rapid spin- 
down phase followed by mainly isotropic evaporation at foreseeable energies. Although the 
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Hawking temperature does not depend strongly on the angular momentum of a spinning 
black hole of a given mass, there is a strong bias in the emission spectra towards higher- 
energy emissions into higher partial waves, which help the black hole to shed its angular 
momentum. The resultant spectra are flatter, with substantial tails beyond 1 TeV. As 
a consequence of these more energetic emissions, rotating black holes emit with reduced 
multiplicity relative to their non-rotating counterparts. However the absolute multiplicity 
can still be large. 

The preferential equatorial emission of scalar, fermionic and high energy vector parti- 
cles leads to slightly less central distributions at detector level. This effect is reduced by 
the evolution of the spin axis during evaporation (away from the initial orientation per- 
pendicular to the beam direction). The emission of polarised higher-spin fields is favoured, 
compared to the spinless case, leading to increased vector emission and marking a fur- 
ther departure from a purely democratic distribution of particle species. This shows little 
dependence upon the number of dimensions. 

These findings will complicate the interpretation of black-hole events, should they oc- 
cur at the LHC or future colliders. While the basic signature of energetic, democratic 
emission of all Standard-Model species and large missing energy remains valid, the deduc- 
tion of the fundamental Planck scale and the number of extra dimensions will be more 
difficult than was anticipated in earlier studies [39]. On the other hand, many interesting 
new and potentially observable features emerge, such as the different angular distributions 
and polarisation of particles of different spins. We intend to investigate possible analysis 
strategies in a future publication. 

In view of the important effects of angular momentum, we would counsel against the 
use of black-hole event generators that neglect these effects. This leaves BlackMax and 
CHARYBDIS2 and the generators of choice for future studies. Both programs take black-hole 
angular momentum fully into account, but they have other features and emphases that are 
complementary. In the formation phase, BlackMax uses a geometrical approximation for the 
cross section and parametrizes the loss of energy and angular momentum as fixed fractions 
of their initial-state values, whereas CHARYBDIS2 incorporates a more detailed model based 
on the Yoshino-Rychkov bounds and comparisons with other approaches. The treatment 
of the evaporation phase in the two programs appears broadly similar, but BlackMax has 
options for brane tension and split branes, and for extra suppression of emissions that 
would spin-up the black hole, while CHARYBDIS2 includes treatment of the polarisation of 
emitted fermions and vector bosons. The conservation of quantum numbers is also treated 
somewhat differently. At the Planck scale, BlackMax emits a final burst of particles with 
the minimal multiplicity needed to conserve quantum numbers, whereas CHARYBDIS2 has 
a wider range of options. 

A deficiency of both programs is the absence of gravitational radiation in the evapora- 
tion phase. This is because the greybody factors have not yet been computed for this case, 
due to extra theoretical difficulties in the separation of variables. Unlike Standard-Model 
particles, gravitons will necessarily be emitted into the bulk, giving rise to a new source 
of lost energy and the possibility of recoil off the brane. In the non-rotating case, it is 
known [73, 74] that bulk graviton emission is small for low numbers of extra dimensions. 
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but increases rapidly in higher dimensions due to the growing number of polarisation states. 
However, the large number of Standard-Model degrees of freedom ensures that brane emis- 
sion remains dominant. Clearly a full treatment of the rotating case is desirable, but there 
is hope that the effects will not be too significant, taking into account the uncertainties in 
energy loss already allowed for in the formation phase. Furthermore, as already discussed 
in Sect. 3, we expect the residual SM charges of the black hole to prevent recoil off the 
brane in all but a small fraction of events. 

In summary, the simulation of black hole production and decay at hadron colliders 
has seen rapid advances in recent years, but there remain substantial challenges in both 
theoretical understanding and data analysis, should such events be seen. It is hoped that 
CHARYBDIS2 will serve as a convenient basis for continuing theoretical refinement and as a 
useful tool in the design of analysis strategies for this exotic and complex possibility. 

Acknowledgements 

We thank colleagues in the Cambridge SUSY Working Group for helpful discussions. 
SRD and MC gratefully acknowledge financial support from the Irish Research Council 
for Science, Engineering and Technology (IRCSET). SRD and MC would like to thank 
Elizabeth Winstanley and Panagiota Kanti for helpful advice and guidance. MC was 
partially funded by Fundagao para a Ciencia e Tecnologia (FCT) - Portugal through 
PTDC/FIS/64175/2006. MS was supported by the FCT grant SFRH/BD/23052/2005. 

Appendices 
A. Conventions 

Throughout this article, unless stated otherwise, all theoretical expressions are in natural 
units where the Planck constant and speed of light are respectively h = c= 1. We keep the 
dependence on the (4 + ra)-dimensional Planck mass explicit in all expressions and adopt 
as reference convention, the PDG definition^^ which uses the Einstein-Hilbert action 

Seh = ImI+- J ^i('+")x^7^(4+„) = J c^(^+")x^7^(4+„) , (a.i) 

to set the reduced Planck mass Md- The Planck mass Md is then defined as 

M|+" = (27r)"M2+" . (A.2) 
An alternative convention is obtained by defining 

M|+;^ = 2M2+" . (A.3) 

This is the Giddings-Thomas convention [6] used internally in CHARYBDIS2. 

'^^See for example the extra dimensions section of the PDG review [105] 
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B. The Yoshino-Rychkov mass/angulcir momentum bounds 

This appendix is divided into two sections. The first is a brief review of the method 

in the Yoshino-Rychkov paper [22], focusing on the key equations we used to implement 
the Yoshino-Rychkov boundary curves in CHARYBDIS2. The second section explains how 
CHARYBDIS2 calculates the boundary curve £,b{C) for a given D and b. 

Summary of the Yoshino-Rychkov method 

It should be noted that in the Yoshino-Rychkov calculation, the usual assumptions are 
made - namely, that the extra dimensions may be regarded as infinite in extent, and that 
the only effect of the branc's gravitational field is to restrict the black hole produced to lie 
on the brane (the 'probe brane' approximation). 

The first stage of the calculation is the choice of metric for each of the colliding partons. 
Yoshino and Rychkov use the Aichelburg-Sexl metric [107]. This is the metric produced 
by boosting a Schwarzschild-Tangherlini metric to the speed of light, whilst reducing its 
associated mass to zero in such a way that its energy remains finite. It is appropriate 
to an ultrarelativistic point particle with no spin or charge, and so the Yoshino-Rychkov 
calculation neglects the effects of spin, charge, and finite size of the partons. 

In the spacetime outside the future lightcone of the collision event (i.e. regions I, 
II and III of Fig. 21), the metric for the complete system is obtained by combining two 
Aichelburg-Sexl metrics corresponding to partons travelling in opposite directions (in the 
centre of mass frame). This gives the correct spacetime outside region IV of Fig. 21, 
because the colliding partons are taken as travelling at the speed of light, so there can be 
no interaction between their gravity waves before the collision. 



Collision 
event 




Left parton and Right parton and 

accompanying shock accompanying shock 

Figure 21: Spacetime regions in a parton-parton collision. In this diagram, the z axis is defined to 
lie along the direction of motion of the left parton, and D — 2 spacelike dimensions are suppressed. 

The next stage in the calculation is the selection of a spacetime slice somewhere in the 
union of regions I, II and III, on which one will look for an apparent horizon (AH). An 

AH is a surface whose outgoing null geodesic congruence has zero expansion. Assuming 
the cosmic censorship hypothesis [108], an event horizon (EH) must be present outside any 
AH; thus finding an AH is sufficient to show that a black hole forms. Furthermore, the 
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fact that the EH must he outside the AH, combined with the area theorem [109] which 
states that the EH area never decreases, can typicaUy be used to set bounds on the mass 
and angular momentum of the formed black hole. 

The slice used by Yoshino and Rychkov is the futuremost shce outside region IV - i.e. 
the boundary of this region. This slice was chosen because it gives the most restrictive, and 
therefore best, bounds on the maximum impact parameter for black hole formation b^ax 5 
and thus on the cross section a = Trb"^^^. It also gives the best bounds on the mass and 
angular momentum trapped in the black hole following production. The reasoning behind 
these statements is detailed more fully below. 

To obtain a bound byRmax on the maximum impact parameter for black hole formation 
for a given D, the impact parameter is increased at the given D until one no longer finds 
an AH - the critical value of b at which an AH no longer appears is then byRmax- This will 
be a lower bound - there may be impact parameters greater than byRmax for which an AH 
forms, but only 'after' the slice considered. 

Prom the above, the best lower bound on bmax that can be achieved (given that we 
are not able to look for AHs in region IV) is obtained by using the futuremost slice outside 
region IV - i.e. the Yoshino-Rychkov slice. The bounds that Yoshino and Rychkov have 
obtained for bmax, and the corresponding results for a, are given in Table 11 of [22]. 

We now discuss the calculation of the mass and angular momentum bound for a given 
b and D in the Yoshino-Rychkov method. In a trapped surface method, this is achieved by 
calculating the D — 2 dimensional area corresponding to the AH, A ah- Since the AH has 
the true black hole EH outside it, and the black hole EH area can never decrease according 
to the area theorem, it is normally the case that the D — 2 dimensional area of the final 
produced black hole EH, Aeh, should be greater than Aah- 

Aah < Aeh ■ (B.l) 

Now, we expect the horizon area of a black hole to be linked to the mass contained within, 
so we can convert the above formula into one in terms of mass. We define the AH mass as 
the mass of a Schwarzschild-Tangherlini black hole with area Aah- 



In (B.2), Gd is the D-dimensional gravitational constant as in Eq. (A.l), and is the 

{D — 2)-area of a unit sphere, given by 

p+i 

With the definition of the AH mass given in (B.2), the equivalent to (B.l) must be 

MAH<Mirr , (B.4) 

where Mj^r is the irreducible mass of the produced Myers-Perry black hole — this is the 
mass of a Schwarzschild-Tangherlini black hole having the same horizon area as the Myers- 
Perry black hole. Since it is defined in terms of the area of the Myers-Perry black hole it 
is a function of both M and J. 
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Equation (B.4) then represents the trapped surface bound on the mass and angular 
momentum trapped in the black hole during production, with the equation of the boundary 
itself being characterised by 

Mah = Mirr . (B.5) 

To convert (B.5) into a boundary line in the (M, J) plane (which can be scaled to a 
boundary line in the (^, () plane, using notation from Sect. 2.2) we first need an equation 
for the irreducible mass of a Myers-Perry black hole with mass M and angular momentum 
J. This may be extracted from the definition of Mirr- 



AMyers-Perry{M, J) — Aschwarzschildi^irr) — ^D-2rg (M, 



(B.6) 



Computing the left hand side of (B.6) using the Myers- Perry metric [47], we find the link 
between M, J and M^rr of a Myers-Perry black hole: 



r^-\Mirr) = r^-HMMM, J) . 
In the above, rniM, J) is the Myers-Perry horizon radius, given by 

'{D-2)Jl'^ 



2M 



r^-'{M)r%-''{M,J), 



5-Dt 



(B.7) 



(B.8) 



whilst rs{M) is the horizon radius of a Schwarzschild-Tangherlini black hole of mass M, 
given by 



rs{M) 



l/(D-3) 



(B.9) 



{D-2)^D-2_ 

We now combine equations (B.5) and (B.7) to show that a point on the bound with mass 
M has horizon radius rnbi^) given by 



VHbiM) 



(Mah) 



„D-3 



(M) 



(B.IO) 



We insert this result into equation (B.8), and then make repeated use of equation (B.9) 
to replace all rs symbols in the result by the explicit expression for this quantity. After 
some rearrangement, we discover that a point on the bound with mass M has angular 
momentum where 



Jb{M) 



2Mah 



IGttGdMah 



{D-2) Y{D -2)^0-2 



l/(D-3) 



{M/Mah 



\D-2 



1 . 



(B.ll) 



This is an explicit equation for the trapped surface boundary line in the (M, J) plane. One 
can rearrange (B.ll) to make M the subject, to obtain an alternative form for the equation 
of the boundary line: 



Mfe(J) = Mah < 1 + 



(L» -2)J' 
2Mah 



1 2 



(D - 2)17g-2 
WttGdMah 



2/(D-3) 



) 



l/(D-2) 



(B.12) 
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To produce valid bounds from the approach of finding an (M, J) bound outlined above, 
it is necessary that an arbitrary surface outside the AH in the slice used has a larger 
area. Unfortunately, for the particular case of the Yoshino-Rychkov slice, this does not 
hold. However, Yoshino and Rychkov found a different area. An,, which they demonstrated 
would be a true lower bound on Aeh (^/6 is equal to twice the area of the intersection of 
the AH with the transverse collision plane). The calculation of the (M, J) bounds then 
proceeds as described above, with Mn, and An, replacing Mah and Aah in equations (B.l) 



The Yoshino-Rychkov slice should give the best bounds in this case as well. This is 
because the use of the futuremost slice outside region IV should give the largest AH areas 
obtainable outside this region [22]. As a result, one expects to get the largest, and therefore 
strictest, lower bounds on Mj^y using this slice. 

Calculation of the mass/angular momentum boundary 

If we rewrite equation (B.12) in terms of the fractions of initial state mass and angular 
momentum ^ and (, replacing all Mah symbols in the equation by Mn, (in accordance with 
the fact that the Yoshino-Rychkov calculation uses the latter quantity), then we obtain 



In (B.13), stands for the energy of each colliding parton in the centre of mass frame. 
This equation can be used to calculate the Yoshino-Rychkov bound ^biC) ^or a given b and 
D, provided one is able to obtain ^n, for the b and D values used. 

The problem of implementing the Yoshino-Rychkov bound in CHARYBDIS2 is then one of 
ensuring that the program has a means of obtaining for all values of D and b (5 < D < 11, 
< b < bYB.max{D)). Calculating ^^f^ exactly (i.e. using the full analytic method) for a given 
b and D, requires repeating the entire Yoshino-Rychkov calculation each time - i.e. find 
the AH numerically, compute the area of its intersection with the transverse collision plane 
All,, and then convert this to an ^n, value using an equation similar to (B.2). Implementing 
the full calculation in CHARYBDIS2, to be repeated on an event by event basis, would slow 
down the program considerably. Instead, we use the ^ib vs. b data files for D = 5 to D = 11 
that had already been generated by Yoshino and Rychkov, and which were used to produce 
the = 5 to D = 11 plots in figure 10 of [22]. CHARYBDIS2 computes the value of for a 
given b and D by linear interpolation on the points in these data files. Linear interpolation 
provides sufficient accuracy due to the close spacing (in b) of the points in the data files. 

C. Details of the 'constant angular velocity' bias 

This section describes the implementation of the constant angular velocity bias in the 
simulation of the black hole production phase. With the bias off (CVBIAS=.FALSE. ), the 
values of (^jC)^*^ are simply those generated from the linear ramp distributions described in 

are the trapped mass and angular momentum fractions passed to the routine which imposes the 
Yoshino-Rychkov boundary condition. 



- (B.12). 




(B.13) 
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Sect. 2.2. When the bias is turned on, the (") point to be passed to the Yoshino-Rychkov 
boundary routine is obtained in a more complex fashion, which is outhned below. 

First, a point is generated using the linear ramp distributions as before. The horizon 
angular velocity and a* value corresponding to the point, Q and a*(^, Q, are calculated 
using the standard equations, and compared to those of the initial state, 0(1, 1) and a* (1, 1). 
In particular, the quantities 10(^,0 — 0(1,1)1/0(1,1) and | log[a*(^, C)/a*(l, 1)]| are cal- 
culated, and compared to the values of the constants A and A respectively (whose values 
win be discussed shortly). Note that |0(^, C) - 0(1, 1)|/0(1, 1) and | log[a*(^, C)/a*(l, 1)]| 
are essentially both measures of the differences between the values at the point and the 
initial state values. 

The point is then assigned a number a(^,C) between and 1, whose value largely 
depends on whether both |0(^, C) - 0(1, 1)|/0(1, 1) < A and | log[a*(^, ()/«*(!, 1)]| < A or 
not (i.e. whether 0(,^, (") and a*(^,(^) are sufficiently close to 0(1,1) and a*(l,l) or not). 
If one or both of the conditions are not satisfied, then the point is assigned a constant 
A; < 1 (as defined below). If both conditions are satisfied, the point is assigned the value 
of a function x(^)C)- The function x(^,C) has the key properties k < x(C,C) ^ !> ^^id 
approaches 1 as 0(^, () gets closer to 0(1, 1). The details of our choices for k and x(^, Q 
will be discussed shortly. 

A random number /? is generated according to a uniform distribution between and 1. 
If a(4, C) > /5 then the point is accepted, otherwise it is rejected. If the point is rejected, 
further (^, () points have to be generated by the ramp distributions, and put through the 
above procedure, until a point is accepted. The final point is passed to the Yoshino-Rychkov 
boundary routine. 

It is reasonably clear that this procedure for generating a (^, () point (to be passed 
to the Yoshino-Rychkov boundary routine) is equivalent to a procedure which generates a 
point from a biased probability distribution of the form asserted in Sect. 2.2. To be specific, 
the biased probability distribution resembles the basic ramp probability distribution, but 
all points whose O and a* values are sufficiently close to those of the initial state have had 
their probabilities enhanced. The enhancement is greater the closer 0(^, (^) is to 0(1, 1). 

We now discuss our choices for the function and the parameters used in the above 
procedure. A suitable choice for the function x(C)C)) which has the properties stated, is 
based on the Breit-Wigner form (note that we introduce a further 'width' parameter T): 

rV4 

^ ([0(^,0- 0(1, i)]/o(i,i))2 + r2/4 ■ ^^-^^ 

The constant k may then be fixed by imposing continuity on a(^,("), such that the 
biased probability distribution represented by the above procedure does not possess any 
sudden jumps. Note that we hope the dividing curves between the enhanced region and the 
unenhanced regions to be |0(^, (^) — 0(1, 1)|/0(1, 1) = A on either side of the 'connected' 
curve 0(^, (^) = 0(1,1) (which is connected to the point ^ = 1,C = 1)- As explained in 
the main text, the a* condition is only present to remove probability enhancement around 
the other 'disconnected' 0(^, (") = 0(1, 1) curve, and should not interfere with the angular 
velocity based enhancement around the right curve. 
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On the basis of these assumptions, continuity of a(^, Q is assured by taking 



k = x(^,C)l|n(«,c)-n(i,i)|/n(i,i)=A 



r2/4 



(C.2) 



A2 + r74 ■ 



Note that with this choice of x(^, Q and fe, F is a constant which sets the 'width' of the 
probabihty peak around the connected Q = 1) curve, and the 'height' of this 
peak. The smaller T is, the sharper and stronger the probability enhancement around the 
appropriate curve. The default value of T is 0.4 for a probability enhancement which is 
not too strong - however, the value of this variable could potentially be changed. 

The values of A and A were chosen by looking at a large number of individual (6, D) 
cases used by CHARYBDIS2, and trying to find a suitable combination of values that gave 
enhancement of a suitable region only around the connected curve. This procedure resulted 
in the choice A = 0.2 and A = 0.4 (these values should not be changed, as small changes 
can cause drastic changes in the regions where the probability is enhanced). 

A final comment is appropriate explaining the slightly peculiar form of the a* condi- 
tion for probability enhancement - | log[a*(^, 1)]| < A. The reason for this form is 
because, when trying to find a condition that discriminated points around the connected 
curve from those around the disconnected curve, we noticed that those around the discon- 
nected curve, had a* values that were one order of magnitude (or more) away from the a* 
values of the points around the connected curve. To remove the enhancement around the 
disconnected curve, a condition based on the logarithm of the ratio a*(^, ^)/a*(l, 1) is then 
more appropriate. 

D. Evaluation of spheroidal wave functions 

In a paper by Leaver [92] the spheroidal wave functions for arbitrary spin are expanded as a 
series inx = cos 9 around x = —1. In this paper there is an extra parameter which can take 
either sign. We can alternatively expand around x = +1. So in general, we can construct 
three more expansions. Even though all four of them converge uniformly in x £ [—1, 1], 
the numerical errors behave differently if we consider a fixed region. Therefore it is useful 
to look at all the options and use different expansions for different regions. The only extra 
complication will be to match them in a common region. 
The possible expansions are 



where a, /3 are chosen as to reproduce the correct behaviour around the regular singular 
points X = ±1, 




(D.l) 



n=0 



m — h 



m + h 



(D.2) 



a = 



2 



2 



and s = ±1 for expansions around x = —1 and x = 1 respectively. Substituting into (3.12) 
we obtain the recurrence relation (for n > 0), 



bn 



\2 n + a n[n + a)J V nJ n + a 



and the simplifying condition = ^ p = ±c. We can set 6_i = and choose the 
normalisation = 1. The parameters above are given by the following expressions: 

cr = a(s + 1) +/3(1 - s) 

a(l - s) + /3(1 + s) - 1 - Asp 



ei 



2 



Ak + h{h + l) + c^ -{a + f3){a + 13 + 1) 
e2 = ^ h a + 

+p [q;(s + 1) + /3(s - 1) - s + s/isign(p)] 
7 = a + /3- l + /isign(p) . (D.4) 

For almost all values of Ak-, the ratio bn+i/bn in (D.3) goes to 1/2 as n — oo. So as the 
order increases, the remainder will behave as (when N — oo) 

J2bn{l + sxr^bN{l + sxf—— , (D.5) 

which diverges at x = s and goes exponentially faster to zero as we move closer to x = —s. 
So if we do not tunc to particular values of Ak this is how the numerical errors will behave. 
However, from the analytical point of view we still want to have an expansion which 
converges uniformly, which is not the case in general as suggested by (D.5). Therefore we 
need to know how the ratio bn+i/bn behaves for large n. Following [110] we determine this 
by assuming that for large n 

bn+i bji 



bn bn-l 



k{n) + ... (D.6) 



where . . . denotes subleading contributions. Inserting this in (D.3) and solving for k we 
get two possible behaviours 

f i -I- 2sp 1 



2sp 
n 







It is straightforward to check that uniform convergence for x G [—1, 1] is only achieved in 
the second case. A solution with these convergence properties is called a minimal sequence 
solution [92]. Furthermore, there is a theorem (see Theorem 1.1 of [111] and the reasoning 
in [110]) which ensures that the sequence obtained from the eigenvalue problem in (D.3), 
with the initial conditions given above, is a minimal sequence. Therefore, after determining 
the eigenvalue wc can compute the solution (up to a normalisation) by fixing 6o = 1- 

However, from the numerical point of view, a small error in the eigenvalues implies that 
bn+i/bn will fail to go to zero when numerically evaluated through (D.3). So the remainder 
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will actually behave as (D.5), since in practice we are approximating the eigenfunction by 
a nearby function, for which the sequence of expansion coefficients behaves as a dominant 
solution of (D.3) when n is large^^. Thus the numerical radius of convergence will be a bit 
smaller and the expansion will fail at a; = s. 

We avoid the above convergence problem by using two expansions, one around each 
singular point x = ±1, and then matching the normalisation in a region where both 
converge appropriately. A simple procedure for matching follows from the observation that 

S+ = AS- 

^log|A| =log|5+|-log|5_| , (D.8) 

where ^4 is a constant and S± denotes expansion around x = ±1 respectively. We can find 
1^1 by averaging the quantity on the right-hand side over an overlap region. Furthermore 
we can estimate the relative matching error through the quadratic deviation, 

^ = Alog|A| 



log^MD-aogl^ll)^ (D.9) 

where (. . .) denotes average over the points used. In practice we use points in the overlap 
region x e [-0.25,0.25]. 

For the expansions above to be useful we need to estimate an order of truncation, 

'T'trunc- A simple Condition is that the ratio of consecutive terms 6„+i/6„ is approximately 
constant above the order of truncation. This is equivalent (in the region sx < where 
there is strong convergence) to an exponentially suppressed upper bound on the remainder 
when ntrunc is large (see (D.5)). From this and (D.3) 

ntrunc ~ max{[7- + \h\ + 4c], [^/\Ak\ + \h\{\h\ + 1) + (c + \m\ + 2\h\ + 3)2]} . (D.IO) 

An initial estimate of the truncation error is obtained from (D.5). Another criterion for 
truncation is that the first neglected higher order term is small compared to the truncated 
sum. 

In practice, after the spheroidal function is calculated with a certain truncation, then 
ten more terms in the series are included and the two estimates are compared. Simultane- 
ously, the first neglected term is compared with the estimate for the sum. If these errors 
are still large then ten more terms are calculated. This procedure is repeated until the 
desired accuracy is obtained. 

To compute the eigenvalues Ak efficiently, the recurrence relation can be put in a 
symmetric tridiagonal form by performing a change of basis, 

an = Xnbn (D.H) 



^'Sco for example [111] for a discussion of numerical issues when generating minimal sequences which are 
solutions of three-term recursion relations. 
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such that the coefficient of the bn-2 term for a certain order n is the same as the coefficient 
of the bn-i term at order n — 1. This is possible if 



n(n + (T) / n{n + a) 



y —sp{n + 1+7) y c(n + 1+7) 

were we have taken (without loss of generality) xq = 1 and the convenient choice p = —sc. 
Furthermore we checked that x„ 7^ Vn > 1, so the transformation is well defined. Thus 
the recurrence relation takes the form 



\/cn{n + a){n + 1 + j)an - n ^— + ^i^ On-i 



+ ^ c{n - l)(n - 1 + a){n + 7)a„_2 = -e20„_i. (D.13) 

This is a tridiagonal symmetric eigenvalue problem which can be solved numerically very 
efficiently, by starting with a truncation order (as estimated above) and checking for pre- 
cision by repeating the calculation with some more corrections. 

E. Decomposition of spheroidal waves into plane waves 

It is well known (sec for example [112]) that in the massless limit, the Dirac fields de- 
composes into two independent fields with helicities h = ±1. If we denote each of those 
2-component spinors by ^h-, then their equations of motion in flat space become 

{dt + ha-d)^h = ^, (E.l) 

where a = (cj^, a^, a^) are the Pauli matrices. Then whichever spatial coordinates x we 
decide to use, the field operator will have the expansion (in terms of positive and negative 
energy solutions) 

^h = Y. -r kAV'?.,A(x)e--^* + S[_,V-;.,A(x)e-^*l , (E.2) 

where a/i^;^,6|j_^ are respectively, the usual fermionic annihilation and creation operators 
for particles and anti-particles; A is an unspecified complete set of quantum numbers; '^h^x 
are spinorial normal modes (a; -|- iha ■ &) i^h,!^ = 0, normalised according to 



d^x 'ili,x^h',X' = 5h,h'S{X, A') ; (E.3) 
5 is a Dirac delta function such that 

5;/(A')<5(A,A') = /(A); (E.4) 

A' 

and the sum sign can represent integration for continuous quantum numbers. Similarly to 
Klein-Gordon theory, it is possible to define a scalar product between 2-spinors ip,x using 
a bilateral derivative: 

(V', X) = ij d\ (i;^dtx - StV'^x) ■ (E.5) 
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Then we find an expression for the operators 



(E.6) 
(E.7) 



can also be expanded in a basis of spinors with a different set of quantum numbers 7, 
which may be associated with a different choice of coordinates. The previous expressions 
wiU then give a Bogohubov transformation between operators in one basis and the other: 

_ i^hu iIjz. ..p~"^7'-1 hi. -I- lib,. ^/=~"'"'^^ ih u ..p"^i'-\ h'. 

',7 



ah,x 



7 

E 



The second term of the first expansion, and the first term in the second one wiU be zero. 
This is because they are responsible for particle creation, which does not occur in this 
change of basis, given that we are in the same Lorentz frame in Minkowski space-time. 
This can be seen explicitly for our case of a transformation between plane wave states and 
spheroidal states. First note that in the Kinnersley basis the plane wave spinors take the 
form [113] 



Xr 



Xr 



e 2 cos -T- cos ■ 



e 2 sm sm ■ 



2 cos sin I + e * 2 sin ^ cos | 



e 2 cos sin 2 + e 2 sm -f - cos 



e 2 cos -f cos ^ — e 2 sm sin | 



(E.IO) 
(E.ll) 



where rip = (0p,0p) defines the orientation of the momentum vector p, (f) = (j) — (f)p 
and we have omitted the e^^'^ dependence. As long as we are looking at a fixed plane 



wave, (f) can be shifted by choosing a new origin, which amounts to setting 



0. The 



upper/lower component of Xp ha-s the interesting property of being invariant under the 
exchange 9p 9. Using the asymptotic form (3.4) for the Cartesian coordinates in terms 
of spheroidal coordinates we obtain 



p • X = cos 6p cos 6\/ (ar)2 + + ar sin 9p sin 9 cos , 



(E.12) 



which is again symmetric under the exchange of ^'s. On the other hand, the spheroidal 
spinors in the Kinnersley basis have the form [29] 



-Rtn.ir)S_.,jJ<^,cos9) 
+4^(r)5i^.„(c,cose) 



(E.13) 



We are seeking for a relation between plane waves and spheroidal waves. This can be 
achieved by writing down the general decomposition 



Xp«'''-'' = E^a(p)x^(x), 



(E.14) 
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where now A = {uj,j,m}. But we know that the upper/lower component of the left hand 
side spinor with hehcity it is invariant under exchange of ^'s, so the cx prefactor must be 
proportional to the upper/lower spheroidal function with argument Op. Furthermore we 
put back the (j)p dependence by shifting to obtain 

^.^.p.x = ^ • 5^,,,. „(c, Op)xi(x) , (E.15) 

where the (/)p dependence in ^ ^. ^(c, $7p) is implicit. An integral relation is obtained 
when multiplying by an appropriate spinor, integrating over x and using the normalisation 
condition (E.3) 

j d3xx^'(x)x^e^P-" = CA • S\^^^ic,np)SH,h'S{u - Up) . (E.16) 
Finally, we can use this in the Bogoliubov transformations to obtain 

4,A J d^P ■ S-h,j,m{c,%)al^p (E.17) 

bix^ J dnp.sij^^{c,np)bl^, (E.18) 

where the prefactor is independent of the angular variables. These expressions show how 
a state of a particle/anti-particle with helicity h decomposes into plane wave states with 
the same helicity and momentum orientations J7p. The probability of a certain orientation 
is given by the square modulus of the spheroidal function with spin weight ^h. 

Since masslcss scalar, spinor and vector perturbations all follow the same master equa- 
tion, these conclusions apply similarly to the remaining cases. 
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